Double hard scattering without double counting

Journal of High Energy Physics, Jun 2017

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

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Double hard scattering without double counting

HJE Double hard scattering without double counting Markus Diehl 0 1 2 5 Jonathan R. Gaunt 0 1 2 3 Kay Sch¨onwald 0 1 2 4 0 Platanenallee 6 , 15738 Zeuthen , Germany 1 De Boelelaan 1081 , 1081 HV Amsterdam , The Netherlands 2 Notkestraße 85 , 22607 Hamburg , Germany 3 Nikhef Theory Group and VU University Amsterdam 4 Deutsches Elektronen-Synchrotron DESY 5 Deutsches Elektronen-Synchroton DESY Double parton scattering in proton-proton collisions includes kinematic regions in which two partons inside a proton originate from the perturbative splitting of a single parton. This leads to a double counting problem between single and double hard scattering. We present a solution to this problem, which allows for the definition of double parton distributions as operator matrix elements in a proton, and which can be used at higher orders in perturbation theory. We show how the evaluation of double hard scattering in this framework can provide a rough estimate for the size of the higher-order contributions to single hard scattering that are affected by double counting. In a numeric study, we identify situations in which these higher-order contributions must be explicitly calculated and included if one wants to attain an accuracy at which double hard scattering becomes relevant, and other situations where such contributions may be neglected. QCD Phenomenology - 1 Introduction 2 3 5 6 7 8 9 3.1 3.2 3.3 4.2 4.3 6.1 6.2 6.3 6.4 8.1 8.2 9.1 9.2 10 Summary A Fourier integrals – i – 4.1 Leading order analysis: collinear factorisation 4.1.1 4.1.2 4.1.3 Squared box graph 2v1 graph Combining contributions Leading order analysis: TMD factorisation Subtraction formalism in momentum space Subtraction terms at higher orders Double DGLAP evolution Resummation of large logarithms Dependence on the cutoff scale Multiscale problems Approximation of the intrinsic distribution Collinear DPDs in momentum space Comparison with other work The approach of Manohar and Waalewijn The approaches of Blok et al. and of Ryskin and Snigirev Quantitative illustrations Simplified analytic estimates Collinear parton luminosities 9.2.1 9.2.2 Setup Numerical results 9.3 Production of two scalars Setting the scene 2.1 Contributions to the cross section: power behaviour and logarithms Short-distance limit of DPDs TMDs Collinear DPDs: splitting contribution Collinear DPDs: all contributions 4 A scheme to regulate DPS and avoid double counting Introduction The precise description of high-energy proton-proton collisions in QCD is imperative for maximising the physics potential of the LHC and of possible future hadron colliders. An important issue in this context is to understand the mechanism of double parton scattering (DPS), in which two pairs of partons undergo a hard scattering in one and the same proton-proton collision. In many situations DPS is suppressed compared with single parton scattering (SPS), but this suppression generically becomes weaker with increasing collision energy. For specific kinematics or specific final states, DPS can become comparable to or even larger than SPS. An overview of recent experimental and theoretical activities in this area can for instance be found in [1, 2]. Consider a DPS process pp → Y1 + Y2 + X, where Y1 and Y2 are observed particles or groups of particles produced in two separate hard scattering processes, whilst X denotes all unobserved particles in the final state. A cross section formula has been put forward long ago in the framework of collinear factorisation, where the transverse momenta q and q2 of Y1 and Y2 are integrated over [3]. A corresponding expression has been given in [4, 5] for transverse momentum dependent (TMD) factorisation, where q1 and q2 are 1 small compared with the hard scales in Y1 and Y2. Inside a proton, the two partons that take part in the two hard scatters can originate from the perturbative splitting of one parton. The relevance of this splitting mechanism for the evolution equations of double parton distributions (DPDs) has been realised long ago [6, 7] and studied more recently in [8–10]. However, it was only noted in [4, 5] that the same mechanism dominates DPDs in the limit of small transverse distance between the two partons, and that the splitting contribution leads to infinities when inserted into the DPS cross section formula. These infinities are closely connected with double counting between DPS and SPS in particular Feynman graphs, a problem that had been pointed out earlier in the context of multi-jet production [11]. Different ways of dealing with this issue have been proposed in [12–15] and [16]. As discussed later in this paper, we find that these proposals have shortcomings either of theoretical or of practical nature. In the present work we present an alternative scheme for computing the cross section in a consistent way, including both DPS and SPS (as well as other contributions). Our scheme allows for a nonperturbative definition of DPDs in terms of operator matrix elements, and it is suitable for pushing the limit of theoretical accuracy to higher orders in αs. Our paper is structured as follows. In section 2 we recall the theoretical framework for describing DPS and specify the theoretical problems mentioned above. The short-distance behaviour of DPDs is discussed in section 3, since it is essential for the scheme we propose in section 4. We first show how this scheme works at leading order (LO) in αs, before giving examples of its application at next-to-leading order (NLO) in section 5. Collinear DPDs evolve according to DGLAP equations, and in section 6 we discuss several consequences of this scale evolution. Our scheme is naturally formulated with DPDs that depend on the transverse distance y between the two partons, but we show in section 7 how one may instead use DPDs depending on the transverse momentum conjugate to y. This allows – 1 – us to compare our results with those of [12, 13] and [14, 15], which we do in section 8. Whilst the focus of the present work is theoretical rather than phenomenological, we give in section 9 some quantitative illustrations of our scheme, obtained with a relatively simple ansatz for the DPDs. We summarise our findings in section 10. Some Fourier integrals required in the main text are given in an appendix. 2 Setting the scene In this section we recall theoretical issues originating from the perturbative splitting mechanism in double parton distributions, namely the appearance of ultraviolet divergences in the naive cross section for double parton scattering, the problem of double counting between DPS and single parton scattering, and the treatment of the so-called 2v1 (two versus one) contributions to DPS.1 We also give some basic definitions and results. Throughout this work we use light-cone coordinates v± = (v0 + v3)/√2 for any four-vector v. write v = (v1, v2) for the transverse components and v = √v2 for the length of the transWe verse vector. Since the perturbative splitting mechanism in DPDs leads to issues in the ultraviolet (UV) region, renormalisation plays a crucial role in our context. We work in dimensional space factors (2π)4 into (2π)D. Bare two-parton distributions are then given by regularisation and extend the definitions of [5] to D = 4 − 2ǫ dimensions, changing phase 2π 2π Z dz1− dz2− ei(x1z1−+x2z2−)p+ p|O2(0, z2) O1(y, z1)|p Fbare(xi, zi, y) = 2p+ dy− with twist-two operators Z 1 Oi (y, zi) = q¯ y − 2 zi 1 W † y − 2 zi Γi W y + 1 2 zi q y + 1 z 2 i zi+=y+=0 for quark distributions, where q denotes the bare field. W (z) is a Wilson line from z to infinity along a prescribed path, which we do not recall here. The Dirac matrices Γi select different quark polarisations. Analogous definitions hold for distributions involving gluons, with quark fields replaced by gluon field strength operators. For ease of writing, we omit colour labels on the operators and distributions throughout this work, bearing in mind that different colour couplings are possible for the four parton fields in (2.1). In the process of deriving factorisation, one finds that the proton matrix element in (2.1) needs to be multiplied by a combination of so-called soft factors, which are vacuum expectation values of products of Wilson lines. More information for the case of single parton distributions can be found in chapter 13.7 of [18] and in the recent overview [19]. A brief account for DPDs is given in section 2.1 of [20], and a more detailed discussion will be provided in [21]. The product Fbare × {soft factors} depends on a parameter ζ which regulates rapidity divergences. A scheme in which a soft factor appears explicitly in the cross section formula was presented in [22]. Since the treatment of soft gluons does not 1We follow here the nomenclature of [17]. The 2v1 contribution is referred to as 4 × 2 in [5], as 3 → 4 in [12, 13] (where four-jet production is considered), and as 1 × 2 in [14, 15]. – 2 – (2.1) (2.2) affect parton splitting in any special way, we will not discuss it further in the present paper. Correspondingly, we will suppress the argument ζ in all DPDs. As a final step one performs UV renormalisation, which we assume to be done in the MS scheme. The DPDs obtained from (2.1) are appropriate for TMD factorisation, with the transverse positions zi being Fourier conjugate to the transverse parton momenta that determine the final state kinematics. These distributions renormalise multiplicatively, with one renormalisation factor for each product of a quark field (or gluon field strength) with the Wilson line at the same transverse position and one factor for each pair of Wilson lines at the same transverse position in the soft factors. Denoting the product of these factors with Z, one obtains the final DPD as F = limǫ→0 Z × Fbare × {soft factors} . It obeys a renormalisation group equation which is a straightforward generalisation of the one for TMDs (given e.g. in [19]). The DPDs needed for collinear factorisation are obtained by setting zi = 0 in Fbare and in the associated soft factors, before renormalisation. Setting zi = 0 introduces ultraviolet divergences in the operators O1 and O2, and in the associated soft factors. The renormalised DPDs are then obtained as F = limǫ→0 Z1 ⊗ Z2 ⊗ Fbare × {soft factors} , where Zi renormalises the operators associated with parton i and where the convolution products are in the momentum fractions xi. In the colour singlet channel, where both operators Oi in (2.1) are colour singlets, the soft factors reduce to unity and one obtains the renormalised twist-two operators that appear in single parton densities. Since the operators associated with partons 1 and 2 renormalise independently (both for the TMD and the collinear case) one may choose different renormalisation scales µ 1 and µ 2 in each of them. This is useful when the two hard subprocesses in double parton scattering have widely different hard scales. In particular, one can then approach the kinematics of the so-called underlying event, with a very hard scattering at scale µ 1 and additional jet production at a much lower scale µ 2 (of course µ 2 needs to remain in the perturbative region for our factorisation approach to be justified). With different scales µ 1, µ 2 in the collinear DPDs, we have a homogeneous evolution equation d d log µ 21 Fa1a2 (x1, x2, y; µ 1, µ 2) = b1 X Z 1−x2 dx′′1 Pa1b1 x x1 x ′ 1 1 x1 , αs(µ 1) Fb1a2 (x′1, x2, y; µ 1, µ 2) (2.3) in µ 1 and its analogue for µ 2. For colour singlet DPDs, the kernels Pa1b1 on the r.h.s. are the usual DGLAP splitting functions for single parton densities. In colour non-singlet channels, both the DPDs and the splitting kernels have an additional dependence on the rapidity parameter ζ [21]. Following the notation of [5], the labels ai, bi in (2.3) denote both the species and the polarisation of the partons. The relevant labels are q, Δq, δq for unpolarised, longitudinally polarised and transversely polarised quarks, likewise q¯, Δq¯, δq¯ for antiquarks, and g, Δg, δg for unpolarised, longitudinally polarised and linearly polarised gluons, respectively. Note that the polarisations of the two partons can be correlated with each other, even in an unpolarised proton. – 3 – To simplify our presentation, we will consider the production pp → V1 + V2 + X of two electroweak gauge bosons Vi = γ∗, Z, W . Our results readily generalise to other processes for which DPS factorisation can be established; in the case of TMD factorisation this requires that the produced particles are colour singlets. We denote the four-momenta of the two bosons by qi, their squared invariant masses by Qi2 = qiµ qiµ and their rapidities by Yi = 21 log(qi+/qi−). We work in the proton-proton centre-of-mass frame, taking the proton with momentum p (p¯) to move in the positive (negative) 3 direction. Furthermore we define xi = r Q2 s i eYi , r Q2 s x¯i = i e−Yi , with s = (p + p¯)2. For the phase space of each gauge boson we have s 2 d4qi = dxi dx¯i d2qi = dYi dQi2 d2qi . 1 2 The “naive” cross section formulae (not taking into account the UV problems discussed below) read dσDPS dx1 dx2 dx¯1 dx¯2 d2q1 d2q2 = 1 X × σˆa1b1 (x′1x¯′1s, µ 12) σˆa2b2 (x′2x¯′2 s, µ 22) Z d2y Fb1b2 (x¯′i, y; µ i) Fa1a2 (x′i, y; µ i) (2.7) for collinear factorisation. The combinatorial factor C is 1 if the observed final states of the hard scatters are different and 2 if they are identical. For simplicity we will consider the case C = 1 throughout this paper, unless mentioned otherwise. As explained in section 2.2.1 of [5], there are further contributions with DPDs that describe the interference of different parton species. They can be discussed in full analogy to the contributions given in (2.6) or (2.7), and we do not treat them explicitly in the present work for ease of notation. The hard scattering cross sections σˆi in (2.6) are for the exclusive final state Vi with transverse momentum qi, which must satisfy qi ≪ Qi. By contrast, σˆi in (2.7) is integrated over qi and inclusive for final states Vi + X. At leading order in αs it contains δ functions that enforce x′i = xi and x¯′i = x¯i. The subtractions for collinear and soft regions in σˆi are different in the two factorisation frameworks, but in both cases they lead to a dependence on the factorisation scale µ i that cancels against the µ i dependence of the DPDs, up to powers in αs beyond the accuracy of the calculation. This happens separately for the two hard scatters (i = 1, 2) and by construction works exactly as in the case of single hard scattering. – 4 – (2.4) As was pointed out in [5], the framework discussed so far suffers from problems in the region of small transverse distances between the two partons in a DPD. The leading behaviour of the collinear distributions F (xi, y) at small y can be computed from the perturbative splitting of one parton into two and gives a behaviour like y−2 up to logarithmic corrections. When inserted in the factorisation formula (2.7) this gives a quadratically divergent integral at small y, which clearly signals an inappropriate treatment of the ultraviolet region. As we will review in section 3.1, the short-distance behaviour of the distributions F (xi, zi, y) is less singular but still leads to logarithmic divergences in the TMD factorisation formula (2.6). Instead of using DPDs depending on transverse positions, one may Fourier transform them to transverse momentum space, integrating 1 2 q = (q1 − q2) . – 5 – 1 (2π)4−4ǫ Z Z d2−2ǫy eiyΔ d2−2ǫz1 d2−2ǫz2 e−i(z1k1+z2k2) for TMDs and for TMDs and collinear distributions. For collinear distributions, this transformation must be made before subtracting UV divergences and setting ǫ → 0: with F (xi, y) ∼ y−2 in D = 4 dimensions, the Fourier integral (2.9) would be logarithmically divergent at y = 0. In D = 4 − 2ǫ dimensions this singularity turns into a pole in 1/ǫ, which requires an additional subtraction as we will review in section 7. Rather than being associated with the individual operators O1 and O2, this subtraction is related to the singularity arising when the transverse distance y between the two operators goes to zero. It leads to an inhomogeneous evolution equation for the DPD F (xi, Δ) in momentum space, which has been extensively discussed in the literature [6–10] for the case Δ = 0. Notice, however, that this extra µ dependence does not cancel in the cross section when (2.7) is rewritten in transverse momentum space. Moreover, the additional UV renormalisation of F (xi, Δ) does not remove all UV divergences at the cross section level. The singularity of F (xi, y) at small y translates into a behaviour F (xi, Δ) ∼ log(µ 2/Δ2) at large Δ (see section 7), which gives a quadratic divergence for the Δ integration in the DPS cross section. It is easy to identify the origin of the UV divergences just discussed. Both in the y and Δ representations, one has integrated over the full range of the integration variable and thus left the region in which the approximations leading to the DPS cross section formulae are valid, namely the region where Δ ≪ Qi or, equivalently, y ≫ 1/Qi (see section 2.1.2 in [5]). Outside this region, the DPS approximations are not only unjustified, but they give divergent integrals in the cross section. This points to another problem, namely that of double counting contributions between SPS and DPS. To see this, let us analyse the graph in figure 1a. Since the transverse boson momenta q 1 and q 2 are approximately back to back (up to effects from the transverse momenta of the incoming gluons) it is convenient to introduce the combination (2.8) (2.9) (2.10) (a) (b) (c) only one proton (indicated by the box). (c) A 2v2 contribution to DPS. Here and in the following it is understood that partons emerging from the oval blobs are approximately collinear to their parent proton. A line for the final state cut across the blobs and the produced vector bosons (wavy lines) is not shown for simplicity. For q ≪ Qi the graph gives a leading contribution to dσ/dq2 if the transverse quark momenta in the loops are all of order Qi. This carries the quark lines far off shell, so that this contribution is naturally associated with SPS, with g + g → V1V2 as the hard subprocess. A leading contribution is also obtained from the region where all transverse quark momenta are much smaller than Qi. This region is naturally described as DPS, with two disconnected hard scattering processes qq¯ → V1 and qq¯ → V2 and double parton distributions with perturbative g → qq¯ splittings, indicated by the boxes in the graph. We denote this as a 1v1 (1 versus 1) contribution to DPS, emphasising its close relation to SPS. A double counting problem for this graph obviously arises if one takes the loop integrals in the SPS cross section over all transverse quark momenta (including the DPS region), and likewise if one integrates the DPDs cross section over the full range of transverse positions, which is equivalent to integrating over all transverse momenta in the quark loops. The problem persists if one integrates the cross section over q. Let us now turn to the graph in figure 1b and consider the cross section integrated over q 1 and q2. In the region of large q ∼ Qi, the quark lines at the top of the graph have transverse momenta of order q and are far off shell. The proper description of this region is in terms of a hard scattering qq¯ + g → V1V2, convoluted with a collinear twisttwo distribution (i.e. an ordinary PDF) at the top and a collinear twist-four distribution at the bottom of the graph. For brevity we refer to this as the “twist-four contribution” henceforth. In the region q ≪ Qi, the g → qq¯ splittings are near collinear and the approximations for DPS are appropriate. We call this the 2v1 contribution to DPS, recalling that there is a qq¯+ g subprocess in the graph. Both small and large q give leading contributions to the integrated cross section, and in a naive calculation adding up the twist-four term and the DPS term has again a double counting problem, as well as divergences in each contribution. The naive DPS cross section has a logarithmic divergence at small y, which is seen by inserting the 1/y2 splitting behaviour of only one DPD in (2.7). In turn, the hard scattering cross section in the twist-four term contains a collinear divergence in the form of an integral behaving like dq2/q2 at q → 0, as we will show in section 4.1.2. – 6 – Clearly, one needs a consistent scheme for computing the overall cross section, without double counting and without divergences in individual terms. An intuitive approach for evaluating DPS is to separate the “perturbative splitting” part of a DPD from its “intrinsic” nonperturbative part.2 This has been pursued independently by Blok et al. [12, 13] and by Ryskin and Snigirev [14, 15]. Taking the intrinsic DPD for each proton, one obtains the 2v2 part of DPS, which does not contain any perturbative splitting and is shown in figure 1c. The splitting part of the DPD is explicitly computed in terms of a single parton distribution function (PDF) and a perturbative kernel. This is multiplied with an intrinsic DPD to compute the 2v1 term. Finally, the product of two splitting DPDs is used to compute the 1v1 contribution in the approach of [14, 15], where an ultraviolet cutoff must be imposed to regulate the quadratic divergence we mentioned earlier. By contrast, the authors of [12, 13] advocate to omit this term and replace it entirely with the SPS contribution to the cross section. We are, however, not able to give a field theoretic definition of the “intrinsic” or “nonperturbative” part of a DPD. The consideration of Feynman graphs in the preceding arguments is instructive, but a satisfactory definition should only appeal to perturbation theory in regions where it is applicable. We regard a nonperturbative definition of DPDs as indispensable for a systematic theory approach, for instance for deriving evolution equations and other general properties. The setup we propose in this work defines DPDs as operator matrix elements as described above, containing both splitting and intrinsic contributions. UV divergences in the DPS cross section are avoided by introducing (smooth or hard) cutoffs in the integrations over transverse distances. The double counting problems are treated within the subtraction formalism used in standard factorisation theorems, described in detail in sections 10.1 and 10.7 of [18] and briefly recalled in section 4 here. The subtraction terms that avoid double counting also remove the above mentioned collinear divergence in the twist-four term. A distinction between “splitting” and “intrinsic” contributions to a DPD will be made in the limit of small transverse distances, where it can be formulated in terms of an operator product expansion (see section 3.3), and when making a model ansatz for DPDs at large distances, which is of course necessary for phenomenology. 2.1 Contributions to the cross section: power behaviour and logarithms In preparation for later sections, we now recall some results for the power behaviour of different contributions to the cross section, referring to section 2.4 of [5] for a derivation. We also recall which logarithms appear in the lowest order 1v1 and 1v2 graphs. As already stated, we take the process pp → V1 + V2 + X as a concrete example, but the discussion readily generalises to other cases. The differential cross section dσ/(d2q1 d2q2) in the region qi ≪ Qi can be computed using TMD factorisation. Here and in the following we write Qi to denote the generic size of Q1 and Q2, and likewise for qi. The transverse momenta qi may be of nonperturbative 2The intrinsic part of a DPD may be studied using quark models [23–30], at least in the valence region, or it may be related to the product of two PDFs if correlations between the two partons are neglected. – 7 – (a) (b) (c) associated with DPS but the loop on the r.h.s. does not. (b) Interference between DPS in the amplitude and SPS in the complex conjugate amplitude. (c) 1v1 graph in the region where the quark loop on the left is collinear whilst the one on the right is hard. This corresponds to the region of SPS/DPS interference, with the boxes indicating DPDs with perturbative splitting for the parton pair in the amplitude. HJEP06(217)83 size Λ or much larger. In the latter case, further simplifications are possible by expressing transverse momentum dependent distributions in terms of collinear ones [21], but we shall not discuss this here. The leading power behaviour of the cross section is dσ When qi goes to zero, it should be replaced by Λ. Three types of mechanisms contribute to the leading behaviour, namely DPS, SPS, and the interference between SPS and DPS. Corresponding graphs are shown in figures 1c, 2a and 2b, respectively. As discussed in the previous section, certain graphs contribute both to DPS and to SPS, depending on the kinematics of their internal lines. The 1v1 graph in figure 1 also has leading regions in which one of the loops is hard and the other is collinear. These regions contribute to the SPS/DPS interference, as shown in figure 2c. The double counting problem thus concerns both SPS, DPS and their interference. Note that the SPS graph in figure 2a contributes to the SPS/DPS interference but not to DPS. Both the amplitude and its conjugate in the 1v1 graph contains a loop integral that ∼ behaves like d2k/k2 in the region Λ, qi ≪ k ≪ Qi. When integrated over the full phase space, each loop thus builds up a so-called DPS logarithm, namely log(Qi/qi) when qi > Λ and log(Qi/Λ) when qi < Λ. Whether these logarithms reside in SPS, DPS or their inter∼ ference depends on how exactly one handles the double counting problem. We come back to this issue in section 4.2. Let us now turn to the cross section integrated over q1 and q2, which can be described using collinear factorisation. The leading power behaviour of the cross section, σ ∼ 1/Qi2, is given by the SPS mechanism alone. Several mechanisms contribute at the power suppressed level σ ∼ Λ2/Qi4. These are 1. DPS, which is suppressed because it can only populate the region qi ≪ Qi rather than the full phase space up to qi ∼ Qi, – 8 – (a) (b) for the other. (b) A graph with twist-three distributions for both protons. 2. the interference between SPS and DPS, which is suppressed for the same reason, 3. hard scattering with a twist-four distribution for one proton and a twist-two distribution for the other. Example graphs are figure 3a, as well as figure 1b with the box removed. figure 3b. 4. hard scattering with twist-three distributions for both protons. An example graph is The rationale for considering such contributions is that — whilst being power suppressed compared with SPS — they may be enhanced by higher parton luminosities at small momentum fractions x, or by coupling constants in the relevant hard scattering subprocesses. Let us emphasise that a complete calculation of the cross section at the level of Λ2/Qi2 corrections would be a formidable task, and it is not even established whether factorisation (in particular the cancellation of Glauber gluons) holds at that level. Notice that in collinear factorisation, the SPS/DPS interference term involves collinear twist three distributions for both protons, because the SPS mechanism forces the two partons in the interfering DPS amplitude to be at same transverse position (see section 2.4.1 in [5]). In this sense, mechanism 2 in the above list may be regarded as a special case of mechanism 4, with a disconnected hard scattering in the amplitude or its conjugate (see figure 2b). A full treatment of contributions with twist-three or twist-four distributions is beyond the scope of this paper. We remark however that twist-n operators contain n or less than n parton fields, and that different operators are related by the equations of motion. For a detailed discussion we refer to [31–33]. Twist-n operators with n parton fields were called “‘quasipartonic” in [34] and involve only the “good” parton fields in the parlance of lightcone quantisation [35]. These are exactly the fields appearing in the definitions of multiparton distributions, so that graphs with a double counting issue between higher-twist hard scattering and DPS (or the SPS/DPS interference) involve only quasipartonic operators. The matrix elements of quasipartonic twist-three operators in an unpolarised target satisfy the important selection rule that the helicities carried by the parton lines must balance. This excludes three-gluon operators since three helicities ±1 cannot add up to zero. For quark-antiquark-gluon operators it forces the quark and antiquark fields to have – 9 – opposite chirality, i.e. one only has the operator combination q¯σ+j q, where the transverse index j is contracted with the polarisation index of the gluon. As for non-quasipartonic twist-three distributions in an unpolarised target, one finds that they are absent in the pure gluon sector [36], whereas the corresponding quark-antiquark distributions are again chiral-odd [37]. Since chiral-odd distributions cannot be generated by gluon ladder graphs, they lack the small x enhancement that is one of the motivations to keep higher twist contributions in the cross section. We will therefore not discuss them further in this work. Note that corresponding selection rules do not hold for TMD correlators, where an imbalance in the helicities of the parton fields can be compensated by orbital angular momentum. Let us finally recall the appearance of DPS logarithms in collinear factorisation. The 2v1 graph (figure 1b) has a behaviour dσ/dq2 ∼ 1/q2 in the region Λ ≪ q ≪ Qi, which gives a log(Qi/Λ) when integrated over the full phase space. Depending on how the double counting between DPS and the twist-four mechanism is resolved, this logarithm can appear in different contributions to the cross section. We will discuss this in section 4.1.2. 3 Short-distance limit of DPDs In this section, we analyse the behaviour of DPDs in the limit where the transverse distance between partons becomes small compared with the scale of nonperturbative interactions. In this region, the splitting of one parton into two becomes dominant. Generalising results in [5] we give expressions in D = 4 − 2ǫ dimensions, which are necessary in intermediate steps when constructing a factorisation formula for the cross section. 3.1 TMDs A useful choice of position variables for describing the parton splitting mechanism is with Fourier conjugate momenta3 1 y± = y ± 2 (z1 − z2) , 1 2 k± = (k1 − k2 ± Δ) , 1 2 Z = (z1 + z2) K = k1 + k2 . (3.1) (3.2) The relation between DPDs in position and momentum space reads F (xi, y±, Z) = 1 (2π)2−2ǫ Z d2−2ǫK d2−2ǫk+ d2−2ǫk− eiZK+i(y+k−−y−k+) F (xi, k±, K) (3.3) in D = 4 − 2ǫ dimensions. As seen in figure 4, one can identify y+ (y−) as the transverse distance between the two partons on the left (right) hand side of the final state cut in the DPD. Correspondingly, the transverse momentum difference between the partons on the left (right) hand side of the cut is k− (k+). The splitting singularities of the DPDs thus occur at y± → 0 or k± → ∞. 3The momentum Δ is called r in [4, 5]. 1 2z2 ments. Here and in the following, the line for the final-state cut of the spectator partons is not shown for simplicity. The perturbative splitting contribution Fspl,pt to transverse-momentum dependent DPDs in momentum space has been calculated at leading order in section 5.2.2 of [5]. Generalising these results to D = 4 − 2ǫ dimensions, we have Fa1a2,spl,pt(xi, k±, K) = ′ kj− kj+ (2µ )2ǫ + . . . , (3.4) where j, j′ are transverse Lorentz indices and fa0 (x1 + x2, K) is an unpolarised singleparton TMD.4 The ellipsis denotes a term that involves a TMD for polarised partons in an unpolarised proton and depends on K but not on k±. In position space we then get Fa1a2,spl,pt(xi, y±, Z) = ′ yj+ j y2−2ǫ yy2−−2ǫ µ 2ǫ Γ2(1 − ǫ) + − π1−2ǫ × fa0 (x1 + x2, Z) 2απs Taj0j→′a1a2 x1 + x2 x1 x1 + x2 , ǫ + . . . (3.5) using the Fourier integral (A.2), where the term denoted by an ellipsis depends on Z but not on y±. It is understood that for transverse quark or linear gluon polarisation, both Fa1a2 and the kernel T carry additional transverse Lorentz indices. fa0 (x1 + x2, Z) is the Fourier transform of fa0 (x1+x2, K). The form (3.4) gives the leading behaviour of the DPD for large k± ≫ Λ, and correspondingly (3.5) gives the leading behaviour for y± ≪ 1/Λ. If one inserts these results into the cross section formula and sets D = 4, logarithmic divergences appear at y+ = 0 and y− = 0. To make them explicit we transform variables to d2Z d2y+ d2y− e−i(q1+q2)Z e−iq(y+−y−) with q defined in (2.10). Performing the angular integration in 1 Z π d2y e±iqy yj yj′ y4 = δjj′ Z dy y J0(yq) + δjj′ − 2qj qj′ q2 Z dy y J2(yq) , (3.6) (3.7) 4Compared with section 5.2 of [5], the kernel T jj′ used here has the opposite order of indices jj′ and includes a colour factor, e.g. TF = 1/2 for the colour singlet distribution 1Fg→qq¯. where y stands for y+ or y−, we see that the integral with J0 is divergent at y = 0. Given the range of validity of (3.5) one should impose y ≪ 1/Λ in (3.7), although the integrals are finite for y → ∞ due to the oscillations of the Bessel functions. The perturbative splitting contribution to DPDs at higher order in αs involves graphs with additional partons radiated into the final state as shown in figure 5a, as well as virtual corrections. It is natural to expect that it will again be singular at y± = 0. A calculation of this contribution is outside the scope of the present work, so that we will limit our analysis of TMD factorisation to perturbative splitting at LO. To compute the DPD cross section, we must also consider the case where only one of the distances y+ or y− is small, whereas the other one remains large. In this case, one has a perturbative splitting only on one side of the final-state cut, as illustrated in figure 5b. We will not discuss the detailed expression of the DPD in this regime, but give its general structure. Setting D = 4 for simplicity, we have Fα1α2,y−→0(xi, y±, Z) = ′ Uαj0→α1α2 (xi) ∗ Dα1α2|α0 (xi, y+, Z) , (3.8) ′ j y− y 2 − were Uα0→α1α2 is a kernel for the splitting α0 → α1α2 in the amplitude (hence its complex conjugate appears in (3.8)). Dα1α2|α0 is the position space version of a transversemomentum dependent twist-three distribution, constructed from the hadronic matrix element p φ0 (Z − y+) φ1 (Z + y+) p 1 − 2 Z = p φ0 φ2 1 2 1 2 1 three parton fields. as shown in figure 5c. We have where φi is a “good” field for parton αi (cf. section 2.1). Distributions Dα0|α1α2 (xi, y−, Z) where α0 belongs to the amplitude and α1, α2 to the complex conjugate amplitude are defined in analogy. In the second step of (3.9) we have used translation invariance and shifted the parton fields to the same position as in the corresponding DPD (see figure 4). The labels αi denote the parton species; it is understood that in (3.8) and the following equations parton helicities are taken fixed on the l.h.s. and must be appropriately summed over on the r.h.s. Note the difference between this notation and the labels ai, which denote parton species and polarisation (none, longitudinal, transverse or linear) and thus refer to a pair of parton legs. The notation with ai is hence not suitable for distributions with If y+ is small, then Dα1α2|α0 (xi, y+, Z) itself can be generated by perturbative splitting, Dα1α2|α0,y+→0(xi, y+, Z) = j y + y2+ Uαj0→α1α2 (xi) fα0 (x1 + x2, Z) . (3.10) Notice that a quark and antiquark produced by perturbative splitting have opposite helicities, so that the corresponding quark-antiquark operator φ2 φ1 in Dqq¯|g must be chiral-even. with perturbative splitting only to the right of the final-state cut. The blob denotes a distribution Dα1α2|α0 . (c) Graph for perturbative splitting in the distribution Dα1α2|α0 . HJEP06(217)83 Inserting (3.10) into (3.8) we obtain Fα1α2,y±→0(xi, y±, Z) ′ y +2 y − = yj+ y2j− Uαj0→α1α2 (xi) Uαj′0→α1α2 (xi) ∗ fα0 (x1 + x2, Z) . Taking appropriate linear combinations of parton helicities, we recover the form of Fa1a2,spl,pt in (3.5) at ǫ = 0. 3.2 Collinear DPDs: splitting contribution We now turn to collinear DPDs, i.e. to the case where z1 = z2 = 0. Let us first consider distributions for two unpolarised or two longitudinally polarised partons, so that the DPDs do not carry any transverse Lorentz indices. The lowest order splitting has been computed in [5]. For 4 − 2ǫ dimensions, one finds the general form Fa1a2,spl,pt(x1, x2, y; µ ) = y2−4ǫ π1−2ǫ µ 2ǫ Γ2(1 − ǫ) fa0 (x1 + x2; µ ) αs(µ ) The kernel for the splitting g → qq¯ reads for instance Pg→qq¯ (u, ǫ) = f u2 + (1 − u)2 − ǫ 2 1 − ǫ with a factor f = 1 for the colour singlet and f = −1/√N 2 In terms of the kernel in (3.4) we have Tgj→j′ qq¯ = δjj′ Pg→qq¯. We recognise in Pg→qq¯(u, 0) the − 1 for the colour octet DPD. usual DGLAP splitting kernel without the terms proportional to δ(1 − u). Going beyond leading order, one can deduce the general form of the perturbative splitting contribution using dimensional analysis and boost invariance. For colour singlet distributions one finds 1Fa1a2,spl,pt(x1, x2, y; µ ) = 1 π1−ǫ y2−2ǫ 1 Z X a0 x1+x2 dv fa0 (v; µ ) v v Va0→a1a2 v v x1 , x2 , αs(µ ), yµ, ǫ (3.11) (3.12) (3.13) (3.14) in D = 4 − 2ǫ dimensions. The convolution integral over v is familiar from factorisation formulae for hard scattering processes. Both f and V on the right-hand side are understood to include all necessary subtractions, so that they are finite at ǫ = 0. The splitting kernel V is a double series Va0→a1a2 (v1, v2, αs, yµ, ǫ ) = n X (yµ )2ǫm Va(0n→,ma)1a2 (v1, v2, ǫ) . s (3.15) The µ (and thus on dimensional grounds the y) dependence of V follows from the fact that the mass parameter of dimensional regularisation appears in graphs only via µ 2ǫαs(µ ); terms with n > m in (3.15) are due to the subtractions of ultraviolet or collinear divergences. At lowest order, the hard splitting graphs are disconnected (with no partons across the final state cut), so that V (1,1)(v1, v2, ǫ) = δ(1 − v1 − v2) V (1)(v1, ǫ). Inserted into (3.14) this gives a form consistent with (3.12). Using that at order αsn the poles of highest order are 1/ǫn−1, we find Va0→a1a2 (v1, v2, αs, yµ, 0) = n X logm(yµ ) Va[0n→,ma]1a2 (v1, v2) . s (3.16) ∞ X α n=1 n m=1 ∞ X α n=1 n−1 m=0 in the physical limit ǫ = 0. For colour nonsinglet DPDs one must regulate rapidity divergences, which complicates the preceding result. Taking e.g. Wilson lines along non-lightlike paths introduces additional vectors and changes the analysis of boost properties of the kernel. We will not pursue this issue here. DPDs with transverse quark or linear gluon polarisation carry transverse Lorentz in′ dices. Their perturbative splitting expressions thus have a tensor structure containing additional factors of yj /y compared with the formulae above. At leading order one readily finds from (3.5) and the appropriate splitting kernels that the factor 1/y2−4ǫ in (3.12) is to be replaced with yj yj /y4−4ǫ times a tensor constructed only from Kronecker deltas. 3.3 Collinear DPDs: all contributions Let us now study the small y behaviour of collinear DPDs in more general terms. We start by writing the relation between unrenormalised DPDs in position and momentum space as Z (2π)D−2F (y) = dD−2Δ F (Δ) + Z dD−2Δ [e−iyΔ − 1] F (Δ) , (3.17) omitting all arguments other than y and Δ. The first term on the r.h.s. is a collinear twist-four distribution, independent of any transverse variable. For small y, the second term is dominated by large Δ, so that one can replace F (Δ) by its approximation for large Δ, following the power counting analysis of section 5.2 in [5]. This leads us to write Fy→0(y) = Fspl,pt(y) + Ftw3,pt(y) + Fint,pt(y) , (3.18) where the three terms on the r.h.s. will be described shortly. In D = 4 dimensions, they respectively go like y−2, y−1 and y0, up to logarithmic corrections. Further terms from the perturbative expansion of F (Δ) give contributions to F (y) that vanish like y or faster. xi and ui denote longitudinal momentum fractions. One may also derive the expansion (3.18) from the operator product expansion, without taking recourse to the transverse momentum representation (3.17). In the definition of collinear DPDs one has a product O2(0, z2) O1(y, z1) of operators with z1 = z2 = 0 but nonzero y. This can be expanded around y = 0 in terms of light-ray operators where all fields are at transverse position 0. These operators have twist 2, 3, 4 for the first, second and third term in (3.18), respectively. The spitting contribution Fspl,pt is given by graphs as in figures 4 and 5a and has already been discussed in the previous subsection. The term Ftw3,pt originates from two types of graphs. The first type involves a single perturbative splitting and a quasipartonic collinear twist-three distribution as shown in figure 5b. The second type has two splittings as in figure 4 and a twist-three distribution with one “good” and one “bad” parton field. Given the helicity constraints discussed in section 2.1, collinear twist-three distributions in an unpolarised proton involve a quark and antiquark with opposite chirality (and possibly an extra gluon). As announced earlier, we discard twist-three terms in the following, since they are expected to become unimportant at small momentum fractions x1, x2. Finally, the term Fint,pt contains contributions without any perturbative splitting; we hence refer to it as the “intrinsic” part of the DPD. It can be written as Fint,pt(x1, x2, y; µ ) = G(x1, x2, x2, x1; µ ) + C(· · ·, y; µ ) ⊗ G(· · ·; µ ) + . . . (3.19) where G is a quasipartonic collinear twist-four distribution and C a perturbative splitting kernel, corresponding to graphs as in figure 6. The convolution ⊗ is in the longitudinal momentum fractions indicated by · · · (cf. figure 6a). The first term in (3.19) corresponds to the first term in (3.17) and is the only contribution that does not involve a hard splitting at all. The ellipsis denotes terms with non-quasipartonic twist-four distributions containing three or two parton fields, together with one or two parton splittings. While having the same power behaviour in y, one may expect that at small x1, x2 these terms become less important than the terms with quasipartonic twist-four distributions, which should roughly grow as fast as the square of two parton densities. We now take a closer look at the second term in (3.19). The kernel C can be determined by computing both sides of (3.19) for a given graph. At order αs only “non-diagonal” interactions, i.e. interactions connecting partons 1 and 2 as in figure 6a and b, contribute to C. The ladder graph in figure 6c is independent of y and thus gives identical contributions to the matrix elements Fint,pt and G. As a consequence it does not contribute to C. At this point we can comment on the scale evolution of the different terms in (3.19). The l.h.s. evolves according to the homogeneous double DGLAP equation for DPDs, which describes “diagonal” interactions, either between the partons with final momentum fraction x1 or between those with final momentum fraction x2. By contrast, the evolution of G(x1, x2, x2, x1; µ ) contains both diagonal and non-diagonal ladder interactions [34]. The non-diagonal interactions in the evolution must thus be cancelled by the µ dependence of the term C ⊗ G. At leading order in αs, this dependence comes only from the coefficient function C, which indeed contains just non-diagonal interactions as just discussed. We finally emphasise that an unambiguous decomposition of F (y) into splitting, intrinsic and twist-three parts is only possible in the limit of small y. If y is of hadronic size, neither the operator product expansion nor the notion of perturbative parton splitting make sense. One may however use the short-distance decomposition (3.18) as a starting point for a model parameterisation of DPDs in the full y range. We describe a simple version of this strategy in section 9. 4 A scheme to regulate DPS and avoid double counting In this section, we present a scheme that regulates the DPS cross section and solves the double counting problem between DPS and SPS, as well as between DPS and the twist-four contribution (figure 1b). Before doing so, we discuss some general considerations that motivate our scheme. bing double parton scattering. The following properties are in our opinion desirable for any theoretical setup descri1. It should permit a field theoretical definition of DPDs, without recourse to perturbation theory. This is the same standard as for the ordinary parton distributions in SPS processes. In particular, it allows one to derive general properties and to investigate these functions using nonperturbative methods, for instance lattice calculations. One may object that so far not even ordinary PDFs can be computed to a precision sufficient for phenomenology. However, important progress has been made in the area of lattice computations, and more can be expected for the future. Furthermore, whereas ordinary PDFs are being extracted with increasing precision from experiment, it is hard to imagine a similar scenario for DPDs, because of their sheer number and because DPS processes are much harder to measure and analyse than most processes from which PDFs are extracted. In such a situation, even semi-quantitative guidance from nonperturbative calculations (such as the relevance of correlations of different types) is highly valuable. As already discussed in section 2, the requirement of a nonperturbative definition prevents us from separating the “perturbative splitting” contribution of a DPD in a controlled way for all distances y (or equivalently for all conjugate momenta Δ). 2. To pave the way for increased theoretical precision, the scheme should permit a formulation at higher orders in perturbation theory. Furthermore, the complexity of the required higher order calculations should be manageable in practice. 3. For collinear factorisation, one would like to use as much as possible existing higherorder results for SPS processes, namely partonic cross sections and splitting functions. This means that the scheme should not modify the collinear subtractions to be made in hard scattering kernels, nor the validity of standard DGLAP evolution for DPDs in the colour singlet channel. 4. For TMD factorisation, it is desirable not to modify Collins-Soper evolution and the handling of rapidity divergences. This again allows one to re-use calculations done for SPS, although rapidity evolution for DPS is necessarily more complicated due to the complexities caused by colour [5, 21]. 5. One would like to keep procedures as similar as possible for collinear and TMD factorisation. This will in particular facilitate the computation of DPS processes at perturbatively large transverse momenta in terms of collinear DPDs [21], adapting the well known procedure for single Drell-Yan production [38]. In principle one can use dimensional regularisation to handle the UV divergences that are induced in the DPS cross section by the perturbative splitting mechanism, as is done with the UV divergences that arise in simpler situations such as single hard scattering. However, contrary to that case, the UV divergences discussed in section 2 arise not at the level of individual DPDs but only when two DPDs are multiplied together and integrated over y. This means that if one treats these divergences in dimensional regularisation, only the product of two DPDs is defined in D = 4 dimensions but not the DPDs separately. This possibility was explored in [16]. However, DPDs and their products remain nonperturbative functions at large y, which according to present knowledge cannot be reduced to simpler quantities in a model independent way. In practice, one therefore needs to model or parameterise them at some starting scale. This is more involved for the product of DPDs than for DPDs themselves, as is the practical implementation of scale evolution. We will come back to this scheme in section 8. Ultraviolet regularisation. We define the regularised DPS cross section by multiplying the integrand in the DPS formula (2.6) for measured transverse momenta with Φ(y+ν) Φ(y−ν) and the integrand in the collinear DPS formula (2.7) with Φ2(yν). The function Φ(u) goes to 1 for u ≫ 1 and to 0 for u → 0, and we can restrict ourselves to the case where 0 ≤ Φ(u) ≤ 1 for all u. More specific requirements are given below. Collinear and transverse-momentum dependent DPDs are defined as specified in section 2, without any modifications. Constructed from operator matrix elements, they contain both splitting and non-splitting contributions. They quantify specific properties of the proton and have a simple physical interpretation, with the same limitations as single parton densities. (We recall that a literal density interpretation of PDFs and TMDs is hindered by the presence of Wilson lines and of ultraviolet renormalisation.) Double counting subtraction. To treat the double counting between DPS and other contributions, we adapt the recursive subtraction formalism of Collins, which we briefly sketch now (details are given in sections 10.1 and 10.7 of [18]). Consider a graph (or sum of graphs) Γ that receives leading contributions from a set of loop momentum regions R. An approximation for Γ is then given by Γ ≈ X CR Γ R with CR Γ = TR Γ − TR CR′ Γ . (4.1) X R′<R In each term one integrates over all loop momenta. The operator TR applies approximations designed to work in momentum region R. Subtraction terms avoid double counting the contributions from smaller regions R′ (regions that are contained in R). In these subtraction terms one applies the approximations designed for R and those designed for the smaller regions. One can show [18] that CR Γ then provides a valid approximation in the region R and in all smaller regions, and P R CR Γ gives a valid approximation of the graph in all relevant regions. All approximations discussed here are valid up to power corrections. In our context, we have graphs in which a set of collinear partons split into partons that can be either collinear (as in DPS) or hard (as in SPS). A slight adaptation of the above formalism is required since we compute DPS using DPDs and a regulating function Φ that depend on transverse distances y rather than transverse momenta. A collinear splitting region R′ then corresponds to large y and the corresponding hard region R to small y, but we keep the ordering of regions R′ < R from momentum space when implementing (4.1). We will show in section 4.3 that our use of subtractions in position space is equivalent to the one in momentum space up to power suppressed effects. The subtraction terms for the DPS region turn out to have a very simple form. They can be obtained by replacing the DPDs in the UV regularised DPS cross section with their appropriate limits for small y± in TMD factorisation and for small y in collinear factorisation. Details will be given in sections 4.1, 4.2 and 5. Criteria 1 and 5 above are obviously satisfied in this scheme. Regarding criterion 2, we note that the higher-order calculations required for the double counting subtraction terms are for the short-distance limit of DPDs, which involve much simpler Feynman graphs than the full scattering process. The introduction of a function Φ in the DPS cross section avoids an explicit modification of the definition of DPDs and thus respects criteria 3 and 4. In particular, the collinear DPDs F (xi, y) in transverse position space follow the homogeneous DGLAP evolution equation (2.3). Since for colour singlet DPDs, the evolution kernels are the familiar DGLAP kernels, the associated scale dependence in the cross section cancels by construction against the one of the hard cross sections computed with the same collinear subtraction as for SPS. In our scheme, we have introduced an additional momentum scale ν to separate DPS from SPS and the twist-four contribution. In practical calculations one may take ν equal to the UV renormalisation scale µ in DPDs, but we find it useful to keep it separate in the general discussion. As a minimal requirement, ν must be of perturbative size, so For the single PDFs fa(x, µ 0), we take the MSTW2008LO distributions [48]. In (9.9) we have multiplied the PDFs by a function of the xi that does not affect the DPD at small xi, but smoothly cuts it off near the kinematic boundary x1 + x2 = 1. The function we use is that given in equation (3.12) of [9], with n = 2. We also take a different y dependence for different parton species. For this we use a simplified version of the model in section 4.1 of [39], taking the width h to be x independent (corresponding to h(x1, x2) of [39] evaluated at x1 = x2 = 10−3) and setting each h with q− indices to be the same as the one with q+. Thus we have with ha1a2 = ha1 + ha2 (9.10) For the splitting piece of the DPD we generalise our ansatz in (9.3), choosing an initialisation scale that goes to b0/y at small y but freezes to a constant value b0/ymax when y exceeds a value ymax that marks the transition between perturbative and nonperturbative behaviour. This ensures that the single PDF and αs in the splitting expression are never evaluated at too low scales. A suitable choice of scale is µ y = b0 , y∗ = y p1 + y2/ym2ax . (9.12) Such a prescription is very similar to the b∗ prescription used in TMD phenomeno Here we take ymax = 0.5 GeV−1, which is one of the values considered in the recent TMD study [50]. Using the same parton dependent Gaussian damping as in (9.13), we have Fa1a2,spl(x1, x2, y; µ y, µ y) 1 = πy2 exp − 4ha1a2 y 2 fa0 (x1 + x2; µ y) αs(µ y) Pa0→a1a2 x1 + x2 x1 x1 + x2 2π . (9.13) The coupling αs(µ y) is determined by 3 flavour running, starting with the MSTW2008LO value αs(µ 0) = 0.68183. The single PDFs are obtained by taking the MSTW2008LO distributions at µ 0 and evolving them according to the DGLAP equations for nf = 3. To obtain the splitting and intrinsic DPDs at scale µ , as in (9.8), the input forms just discussed must be evolved, starting from µ 0 for Fint and from µ y for Fspl, according to the homogeneous double DGLAP equations. For this we use a modified version of the code developed in [9]. The modified code works on a grid in the xi, µ and µ y directions (the grid of the original code is in xi and µ only). The grid points in the xi directions are evenly spaced in log(xi/(1 − xi)), whilst those in the µ and µ y directions are evenly spaced in log µ or log µ y. The integrals appearing in the double DGLAP equations are computed from points in the xi grids using Newton-Cotes rules (for details see appendix A of [9]), and evolution from one point in the µ grid to the next is carried out using the Runge-Kutta method. – 56 – with 60 points in the µ direction, spanning µ 0 < µ < 170 GeV, and with 60 points in the µ y direction, spanning b0/ymax < µ y < 340 GeV. According to the studies made in [9], this suggests an error on the level of a few per cent in the DPD values obtained after evolution, which is tolerable in this first study. The DPDs computed on this grid are used together with an interpolation code to produce numerical values for the investigations below. production of a W boson pair). Taking the collider energy to be √ We begin with a study where the scale µ is equal to Q1 = Q2 = 80 GeV (as in the s = 14 TeV, we set x1 and x¯1 in (9.1) to correspond to central production of the first system and x2 and x¯2 to correspond to the production of the second system with rapidity Y (all rapidities refer to the pp centre of mass). This gives x1 = x¯1 = 5.7 × 10−3 , x2 = 5.7 × 10−3 exp(Y ) , x¯2 = 5.7 × 10−3 exp(−Y ) . (9.14) a1a2 b1b2 = uu¯u¯u + u¯uuu¯, a1a2 b1b2 = gggg and a1a2 b1b2 = ud¯d¯u + d¯uud¯. In figure 12 we plot La1a2b1b2 (Y ) in the range 0 ≤ Y ≤ 4 for the parton combinations The first parton combination appears e.g. in ZZ production, the second is important in four-jet production, and the last appears in W +W +. For ease of language, we will refer to these parton combinations as the uu¯, gg and ud¯ channels, respectively. We split the overall luminosity into contributions from 1v1 (Fspl × Fspl), from 2v1 (Fspl × Fint + Fint × Fspl) and from 2v2 (Fint × Fint). We vary ν by a factor of 2 up and down around a central value of 80 GeV, in order to see how DPS alone is affected by variation of this cutoff. For each contribution, the line in the plots denotes the luminosity with ν = 80 GeV, whilst the band is generated from the envelope of the functions with ν = 40 GeV, ν = 80 GeV and ν = 160 GeV. In the case of the 2v2 contribution, the ν scale variation is negligible (as expected from basic considerations, see (9.4)), so this appears as a dashed line in each plot. For the 1v1 contribution of the gg and uu¯ channels, we also plot a band generated by varying the prediction with ν = 80 GeV by a factor 4 up and down. This corresponds to a strictly quadratic cutoff dependence of L, i.e. to the variation of ν by a factor 2 in the naive formula (9.6), where DPD evolution is neglected. Any discrepancy between this band and the actual 1v1 band is therefore due to evolution effects. We do not plot this band for the ud¯ channel: there is no LO splitting process giving ud¯, so that the scale variation (and central value) from (9.6) is zero in this case. We immediately notice in figures 12a and b that the 1v1 contribution is generally much larger than the 2v2 and 2v1 contributions, and that it has an enormous ν variation. To obtain a sensible prediction in these cases, one must include the SPS corrections up to the order that includes the double box graph in figure 1a, so that the associated subtraction term can approximately cancel the ν dependence of DPS. We also notice that at central rapidities in the uu¯ channel, the 1v1 ν variation band essentially fills the band corresponding to a quadratic ν dependence, indicating small evolution effects in this case. By contrast, the 1v1 ν variation band for gg is clearly smaller than the naive expectation, which indicates significant evolution effects. We shall explore the reasons for these differences below. 109 2 108 V e G 2 106 V e G / L 105 ud¯ 2 eVG1011 1012 L 1010 109 gg 2v2 1v1 2v1 1v1× 14 ...4 2 Y (b) 0 2v2 1v1 2v1 with Q1 = Q2 = 80 GeV at √ s = 14 TeV, one with central rapidity and the other with rapidity Y . The parton combinations a1a2 b1b2 are uu¯u¯u + u¯uuu¯ (a), gggg (b) and ud¯d¯u + d¯uud¯ (c). In figure 12c, the 1v1 contribution is small compared to the 2v2 piece and has a small ν dependence. This is because, as already mentioned, there is no LO splitting directly giving ud¯ (generation of a ud¯ pair requires at least two steps, such as u → u + g → u + d + d¯). In this case, there is less of a need to compute σSPS and the subtraction term σ1v1,pt up to the order that contains the lowest-order DPS-type loop (in both amplitude and conjugate). This is fortunate, since this order is two powers of αs higher than for graph 1a (two-step rather than one-step splittings are required in both protons), and the corresponding SPS calculation will not be available for some time. In the gg channel, the 2v1 and 2v2 contributions are of similar magnitude and shape. This is consistent with what has been found in previous phenomenological studies [13, 51] at similar scales and x values. In the uu¯ channel, the magnitude of 2v1 and 2v2 is broadly the same but their shape is somewhat different. In the ud¯ channel, the shape of the two contributions is similar, but 2v1 lies well below 2v2 (and is close to 1v1). We can examine the degree to which the ν variation in the gg and uu¯ channels is reduced by adding the remaining fixed-order terms in the cross section, even without considering a specific final state or having to compute σSPS. This is because only σDPS and the subtraction terms depend on ν. Just as we did for σDPS, we can factor the two partonlevel cross sections out of σ1v1,pt and introduce a “subtraction luminosity” L1v1,pt. This is defined as in (9.1) except that the DPDs are replaced by the fixed-order splitting expression evaluated at scale µ , Fa1a2,spl,fo(x1, x2, y; µ ) = 1 fa0(x1 + x2; µ ) αs(µ ) πy2 For a given parton combination a1a2 b1b2, one can directly subtract L1v1,pt from the DPS luminosity LDPS and compare the ν dependence of the result to that of LDPS alone. Up to multiplication by the appropriate parton-level cross sections, the ν dependence of LDPS − L1v1,pt is equal to that of the full cross section. On the other hand, the overall magnitude of LDPS − L1v1,pt (which can be positive or negative) has no direct significance, since one needs to add the SPS contribution to obtain the physical cross section. We plot LDPS and LDPS − L1v1,pt for the uu¯ and gg channels in figure 13. The lines are the values for ν = 80 GeV, whilst the bands show the variation when ν is varied by a factor of 2 up and down. Note that for LDPS − L1v1,pt, the curve with ν = 80 GeV is always at the very top of the scale variation band. We see that the ν scale variation of LDPS − L1v1,pt is indeed reduced compared to LDPS, with the reduction being much stronger for uu¯ than for gg. The latter is consistent with our previous observation that evolution effects are weaker in the uu¯ channel at central rapidities: if evolution effects are weak, Fspl and Fspl,fo have a similar y dependence, so that (Fspl × Fspl − Fspl,fo × Fspl,fo) is flat in y space and LDPS − L1v1,pt varies weakly with ν. In section 6.2 we argued on generic grounds that the ν variation of the 1v1 contribution to σDPS should be of the same order as the subtraction term σ1v1,pt. This is confirmed in figure 13, where the size of L1v1,pt (evaluated at central ν) can be read off from the distance between the central curve for LDPS and the top of the band for LDPS − L1v1,pt. This finding allows us to sharpen the argument we already made in the discussion of figure 12: not only traction luminosities, shown for the uu¯ (a) and gg (b) channels. Kinematic conditions are as in figure 12. In the case of LDPS − L1v1,pt, the value for ν = 80 GeV is always at the top of the scale variation band. does a large ν variation of the 1v1 contribution to σDPS indicate the need to include SPS and the associated subtraction term in the cross section, but the size of the ν variation may even serve as a rough estimate for the size of these terms. To better understand the behaviour we have seen in the 1v1 luminosities for the uu¯ and gg channels, we now take a closer look at how evolution affects the y dependence of splitting DPDs. In figure 14 we plot Fspl(y; µ ) against y, setting x1 = x2 = 5.7 × 10−3 so that the uu¯ and gg DPDs are the ones used for the point Y = 0 in figures 12 and 13. We also plot the ug DPD, which mixes with uu¯ and gg under evolution. For comparison we show Fspl,fo(y; µ ), as well as the initial condition Fspl,ini(y) = Fspl(y; µ y) given in (9.13), from which Fspl(y; µ ) is obtained by evolution. The distributions are plotted in the range b0/(160 GeV) < y < b0/(40 GeV), i.e. in the range over which the lower integration limit in y is varied when we vary ν between 40 GeV to 160 GeV. In this region, the exponential damping factor in our DPD model is irrelevant, so that Fspl,ini and Fspl,fo only differ by the scales taken in αs and fa0 (x1 + x2). For the initial conditions Fspl,ini we note that the uu¯ and ug distributions are of similar size; the former is initialised by a larger PDF (fg instead of fu) but has a smaller splitting coefficient P (1/2) as we noted before (9.6). By contrast, the gg distribution is much bigger; here both the initialising PDF and the splitting coefficient are large. By construction, all three curves in each plot are equal at y = b0/(80 GeV), when µ y = µ . In all plots, Fspl,ini has a more shallow y behaviour than Fspl,fo. This difference is mainly driven by αs, as the PDFs do not vary so much between µ y and µ at momentum fraction x1 + x2 = 1.14 × 10−2. One observes that the difference between Fspl and Fspl,ini is more significant for gg and ug than for uu¯, i.e. that DPD evolution has a much stronger effect on the former two channels. This is to be expected at small x, since the 1/v behaviour The parton combinations shown are uu¯ (a), gg (b) and ug (c). The vertical grey lines correspond to y = b0/(160 GeV), y = b0/(80 GeV) and y = b0/(40 GeV). of the splitting kernels Pgg(v) and Pgq(v) at small v favours the radiation of a gluon with much smaller momentum than its parent parton, whereas the kernels Pqq(v) and Pqg(v) giving a quark stay finite for v → 0. An interesting point to note is that, whilst for gg and ug the curves for Fspl are more shallow than for Fspl,ini, in the case of uu¯ the Fspl curve is actually steeper. The latter is surprising since — based on the experience with single PDFs at small x — one may expect that forward evolution for y > b0/(80 GeV) would always increase a DPD, and backward evolution for y < b0/(80 GeV) would always decrease it. This indeed happens in the gg and ug channels, whilst forward evolution results in a decrease of the uu¯ DPD. The reason for this is that Fug,spl and Fgu¯,spl, which directly feed into Fuu¯,spl during evolution, are comparatively small. In the case of PDFs, fg is much larger than fu and hence can (c) drive its small-x evolution although the splitting function Pqg does not favour the radiation of low-momentum quarks. The evolution of the gg and ug DPDs is driven by the large distribution Fgg,spl and enhanced by the Pgg splitting function. Note that in all three channels, Fspl has a smaller slope in y than Fspl,fo. This implies that the y integrand for the computation of LDPS − L1v1,pt is positive for y > b0/(80 GeV) and negative for y < b0/(80 GeV). Therefore, the LDPS − L1v1,pt curves for ν = 40 GeV and ν = 160 GeV lie below that for ν = 80 GeV. The curve for the central ν value thus lies at the top of the ν variation band, as we already observed in figure 13. For uu¯, the evolution from µ y to µ turns out to have much the same quantitative effect as adjusting the scale in the fixed order expression from µ y to µ , such that Fspl and Fspl,fo end up extremely close together. This explains why at central rapidities the ν variation of the 1v1 uu¯ contribution in figure 12a coincides almost exactly with the naive prediction from (9.6), which assumes that Fspl depends on y like y−2 (as does Fspl,fo). Furthermore, if Fspl and Fspl,fo are very close to each other, then LDPS − L1v1,pt is much smaller than LDPS, which we indeed see in figure 13a. In the gg channel, the modification of the y slope by evolution is significant. In the range of figure 14, the y dependence is changed from y−2 at the starting scale to approximately y−1.4. One may expect this flattening effect to become even stronger as the x values in the DPDs decrease. The dependence of LDPS on the scale ν should then decrease. This was already anticipated in [14, 15], where studies in the double leading logarithm approximation were performed. Here we investigate this effect using our full LO DGLAP set-up. We consider the gg channel with all x values set equal and study the DPS luminosity as a function of the common x value. We fix the collider energy √s, so that the scale µ = Q1 = Q2 = x s also √ to 7 GeV < µ < varies with x. We make one plot with √ 2800 GeV) and another one with higher collider energy √ s = 100 TeV s = 14 TeV and 5 × 10−4 < x < 0.2 (corresponding and 10−4 < x < 0.02 (corresponding to 10 GeV < µ < 2000 GeV). In our numerical code this requires a new DPD grid with larger µ and µ y ranges, which was generated using 60 points in the µ and µ y directions (as in the original grid). The resulting 1v1 luminosity is shown in figure 15, with bands for the ν variation and its naive version as in figure 12. For comparison we also show the 2v1 and 2v2 luminosities. We see that the ν variation does indeed become progressively smaller compared to the central value as x (and µ ) decreases. At the lowest x values it is much smaller than the naive expectation. At larger x, where evolution has a smaller effect on the DPDs, the ν variation becomes larger: for √ s = 14 TeV the actual and naive ν bands essentially coincide dominate over 1v1. For given µ , this effect is more pronounced for √ towards the right of the plot, where values of the order of x ∼ 0.1 are reached. Towards the left of each plot, where x and µ become small, the 2v2 and 2v1 contributions begin to s = 100 TeV than for s = 14 TeV, since the x values in the former case are smaller. We note that in the uu¯ channel (not shown here) the actual and naive ν bands remain very similar throughout the kinematics of figure 15. For the small x values where the ν variation is small in figure 15, the gg splitting DPD is significantly shallower in y than the fixed-order y−2 form. For √ s = 14 TeV we find that 1015 2 1010 V e G / 105 L 100 10−5 102210 1020 1018 2 1016 eVG1014 /1012 L 1010 108 10160−4 10 √s = 14TeV 10−3 √s = 100TeV 10−2 x (a) μ / GeV 100 10−3 x (b) s is fixed to 14 TeV (a) or to 100 TeV (b). The scales of the DPDs, set to µ = x√s, are given at the top of the plot. a behaviour like y−0.5 is reached around 5 × 10−4. For √ around x = 2×10−3 it reaches y−1, implying a logarithmic dependence on ν in LDP S, whilst s = 100 TeV the y behaviour is like y−1 around x = 4 × 10−4 and like y−0.5 around x = 10−4. In such a regime, the bulk of the contribution to the 1v1 part of LDPS comes from large distances y ≫ 1/ν, where the DPS approximations used to derive σDPS are valid. Also, the 1v1 part of σDPS becomes clearly larger than its ν variation. As we argued earlier, the ν variation should be of the same size as the subtraction term. This is confirmed in figure 16, which is analogous to figure 13b but for the kinematics of the point with µ = 10 GeV in figure 15a. Using our argument 1000 2v2 2v1 1v1× 141..v.14 10−1 1000 2v2 2v1 1v1× 141..v.14 10−2 HJEP06(217)83 7.1 × 10−4, corresponding to Q1 = Q2 = 10 GeV at √ the L axis. s = 14 TeV. Notice the suppressed zero on that σ1v1,pt is of the same magnitude as the order of σSPS that contains the lowest-order DPS-type loop, one may in this situation justifiably make predictions that include the DPS piece but omit the order of SPS just specified, as well as the associated subtraction term. This is encouraging, for instance in the context of four-jet production, where the computation of the relevant SPS order (namely NNLO) is well beyond the current state of the art. Note that lower orders in SPS should be computed and included, if possible. Notice that when the y behaviour of Fspl becomes flatter than y−1, the dominant y than the fixed-order form y−2. region in the 1v1 part of LDPS shifts to values y ∼ 1/Λ, where one cannot compute the splitting DPD in perturbation theory and must rely on a model. Likewise, the 2v1 part of LDPS becomes increasingly sensitive to the region y ∼ 1/Λ as soon as Fspl becomes flatter Another kinematic regime where the 1v1 contribution to LDPS becomes large compared with its ν variation and compared with L1v1,pt is when the two hard systems have a large separation in rapidity. This reduction of the ν variation can be seen as Y approaches 4 in figure 12, but only in the uu¯ and not in the gg channel. In both channels, the effect becomes more pronounced once the rapidity separation of the hard systems is increased beyond 4. To illustrate this, we make plots similar to figure 12 but now with one hard system at rapidity Y and the other at rapidity −Y (rather than one at Y and the other at 0), such that a given value of Y corresponds to a rapidity separation of 2Y . The results are shown in figure 17, where we see that the ν variation in the 1v1 contribution is strongly reduced towards the right hand side of the plots, becoming hardly visible in the uu¯ channel. Also notable is the fact that for large Y , the 2v2 contribution in this channel exceeds 1v1, which strongly decreases between Y = 0 and 2. This reduction in ν dependence can be explained by the fact that at large Y we probe splitting DPDs with one large x and one small x parton. From the point of view of small x logarithms, it is preferable to generate such a configuration by having the 1 → 2 splitting 2v2 1v1 2v1 1v1×14...4 Y (a) V V e e G G 1012 2 1011 /L1010 109 108 gg 2v2 1v1 2v1 1v1×14...4 Y (b) s = 14 TeV, one with rapidity Y and the other with rapidity −Y . The parton combinations a1a2 b1b2 are uu¯u¯u + u¯uuu¯ (a) and gggg (b). at large x, generating directly the large x parton plus a gluon with smaller x, the latter of which splits in a number of stages into smaller x gluons, eventually yielding the small x parton. This increasingly happens with increasing y, since the “evolution distance” between the initial and final scales, µ y and µ , is increased. Thus, evolution again flattens the DPD compared to the naive y−2 expectation and reduces the ν dependence of the DPS luminosities. The effect is particularly drastic in the uu¯ channel, because the lowest-order splitting g → uu¯ is inefficient at generating a pair with very different x values (as Pg→qq¯(v) goes to a constant at small v). Therefore, the repeated splitting mechanism described in the previous paragraph is strongly preferred, even though its last step is penalised by the lack of small-v enhancement in Pgq(v). In figure 18 we plot the uu¯ DPD in the x1 ≪ x2 configuration relevant for Y = 4 in figure 17a. It has a similar y dependence as the ud¯ DPD in equal kinematics, with the curve for ud¯ crossing zero exactly at y = b0/µ (by construction), and the curve for uu¯ crossing zero rather close to this point. In the former case, the lowestorder 1 → 2 splitting process is forbidden, whereas in the latter it is numerically almost irrelevant. The situation for x1 ≫ x2 is the same. To end this section, we study polarised distributions and luminosities. We limit ourselves to the splitting part Fspl of the DPD and to the 1v1 contribution to LDPS. In fact, we have little guidance for modelling the intrinsic part Fint, where a product ansatz as in (9.9) makes no sense. We note that in [39], different ansa¨tze were made for polarised DPDs, and the effects of evolution on these distributions were investigated. It was found that, in many cases, evolution to higher scales leads to a suppression of the polarised DPD with respect to its unpolarised counterpart. which is needed for the point Y = 4 in figure 17a. We initialise the polarised Fspl at µ y using the expression in (9.13) with the unpolarised 1 → 2 splitting functions replaced by their polarised counterparts (given in appendix B of [39]). These distributions are then evolved to the scale µ using the appropriate polarised double DGLAP equations (the required polarised splitting functions are collected in appendix A of [52]). The settings we used for the evolution code are the same as specified in section 9.2.1. Q1 = Q2 = 80 GeV at √ We consider the same scenario as in figure 12, i.e. the production of two systems with s = 14 TeV, one with rapidity 0 and the other with rapidity Y , and now include polarised 1v1 contributions. In figure 19a, we reproduce figure 12a but include 1v1 luminosities for the polarised parton combinations a1a2 b1b2 = ΔuΔu¯Δu¯Δu + Δu¯ΔuΔuΔu¯ and δuδu¯ δu¯ δu + δu¯ δuδuδu¯. For brevity we refer to these combinations as ΔuΔu¯ and δuδu¯ in the following, recalling that Δq and δq indicate longitudinal and transverse quark polarisation, respectively. These polarised luminosities appear in the DPS cross section for double Z production, see section 3.1 of [53] (combinations like uu¯ Δu¯Δu also appear but are not shown here). To avoid cluttering the plots, we omit the unpolarised 2v1 luminosity, reminding the reader that it is of the same magnitude as 2v2. In figure 19b we repeat the exercise for the pure gluon channel, including the 1v1 luminosities for the polarised combinations a1a2 b1b2 = ΔgΔgΔgΔg and δg δg δg δg, referred to as ΔgΔg and δg δg, where Δg and δg respectively denote longitudinal and linear gluon polarisation. We see that the 1v1 luminosities for ΔuΔu¯ and δuδu¯ essentially coincide with the one for uu¯ at central rapidities, with ν variation bands of very similar size. These observations can be explained by the fact that at Y = 0 we initialise all three DPDs using the same expression up to a minus sign, with Pg→qq¯ = |Pg→ΔqΔq¯| = |Pg→δq δq¯| at v = 1/2, and that evolution has a rather weak effect on all three DPDs for central rapidities. This close agreement extends out to higher values of Y for the uu¯ and ΔuΔu¯ curves, owing to the fact that Pg→qq¯ = |Pg→ΔqΔq¯| for general v. By contrast, the δuδu¯ curve clearly falls below 109 2 107 / L105 104 103 102 2v2 uu¯ 1v1 uu¯ 1v1 ΔuΔu¯ 1v1 δuδu¯ Y (a) 1012 1011 1010 2 109 eVG108 /L107 106 105 104 103 2v2 gg 1v1 gg 1v1 ΔgΔg 1v1 δgδg Y (b) (and the unpolarised 2v1 band omitted). the other two for larger Y , because Pg→qq¯ > |Pg→δq δq¯| for v 6= 1/2. In figure 19b, the δg δg luminosity lies well below the one for ΔgΔg, which in turn lies somewhat below the one for gg. This is mostly driven by the differences in initialisation expressions for the DPDs: at v = 1/2 (relevant at Y = 0) the relevant splitting function satisfy for instance Pg→gg : Pg→ΔgΔg : Pg→δg δg = 9 : 7 : 1. Ignoring evolution effects, one would then expect the gg luminosity at Y = 0 to be roughly 1.5 times bigger than the one for ΔgΔg, which in turn should be roughly 50 times bigger than the one for δg δg. This expectation is quite close to the actual luminosity ratios, which are ∼ 4 and ∼ 60, respectively. The remaining difference is due to evolution effects, which increase gg rather considerably, increase ΔgΔg to a lesser extent, and hardly change δgδg. Overall, we see that the polarised 1v1 luminosities can be of the same order of magnitude as the unpolarised ones, so that one must in general take into account all possible polarisation combinations in the DPS contribution, together with the SPS term and subtraction (where the subtraction will also contain both unpolarised and polarised contributions). 9.3 Production of two scalars In this section, we study an explicit example to test our general argument that σ1v1,pt is of similar size as the corresponding perturbative order of the SPS cross section. The process we investigate is an artificial one, chosen for ease of computation. We consider the hypothetical production of two identical massive scalar bosons φ, which couple to quarks with a Yukawa term cφq¯q, where c is a constant. For simplicity we take only one light quark flavour; including further light flavours just multiplies all following results (SPS and subtraction) by nf . We will not compare the subtraction term to the full SPS contribution, but rather to the gg initiated part of σSPS, which is the piece that contains HJEP06(217)83 hence be meaningfully considered by itself. In keeping with the notation of the paper, let us denote the mass of each produced boson φ by Q. Let Y be the rapidity of the diboson system in the pp centre-of-mass frame, sˆ its squared invariant mass, and β = p1 − 4Q2/sˆ the velocity of one boson in the diboson centre-of-mass frame. The gg → φφ part of the SPS cross section can be written as dY dβ dσSPS = xx¯ g(x) g(x¯) = xx¯ g(x) g(x¯) 2β 1 128π 1 − β2 σˆgg→φφ 1 − β2 Z qm2ax Q4 0 dq2 p1 − q2/qm2ax (9.16) HJEP06(217)83 where qmax = βQ p1 − β2 and the gg → φφ matrix element squared, |Agg→φφ|2, is averaged over the spin and colour of the incoming gluons. The momentum fractions x and x¯ in the gluon distributions are obtained from the constraints sˆ = xx¯s and Y = 12 log(x/x¯). The matrix element squared for gg → φφ can be obtained from the gg → HH matrix element squared given in [54] by making the replacement παW mq2 M W2 → c2 . In particular, |Agg→φφ|2 in (9.16) is related to the terms “gauge1” and “gauge2” in [54] according to Agg→φφ 2 = c4αs2 210π2mq4 |gauge1|2 + |gauge2|2 . The right hand side of (9.18) is given in [54] for general quark mass mq. We evaluate this expression numerically for very small mq, using analytic mq → 0 approximations where this is necessary to avoid numerical instabilities. Now we turn to the subtraction term. The two elementary qq¯ → φ cross sections in this term contain the kinematic constraint 2πδ(xix¯i s − Q2), which fixes the momentum fractions xi and x¯i entering the two hard subprocesses at given x, x¯ and β up to a two-fold ambiguity. The DPS subtraction term is given by dσ1v1,pt = dY dβ π2 1−β2 Q8 x x¯ 2 2 16 d2y Φ2(yν) X R=1,8 1 RFa1a2,spl,pt 2 x(1+β), 2 x(1−β), y RFb1b2,spl,pt 2 1 x¯ (1−β), x¯ (1+β), y 1 2 +{β → −β} +{F → I} , (9.17) (9.18) (9.19) where we sum over all possible colour, spin, and quark number interference/correlation possibilities. The index R on F denotes the colour channel (R = 1 for colour singlet, R = 8 for colour octet), and in the sum over a1a2 b1b2 we sum over both unpolarised (e.g. qq¯q¯q) and polarised (e.g. ΔqΔq¯Δq¯Δq) combinations. By {F → I} we denote the same expression, but with the quark number diagonal distributions replaced by the quark number interference ones (see section 2 of [5]). The splitting kernels for the relevant spin combinations are [5] 1 2 1 Pg→qq¯ (v) = − 1Pg→ΔqΔq¯ (v) = v2 + (1 − v)2 , 1 Pg→δqδq¯ (v) = − δjj′ v(1 − v) (9.20) and 8P = − 1P/pNc2 − 1, where j, j′ are indices for transverse quark polarisation. The term with interference DPDs I gives the same contribution as that with F , since the corresponding diagrams are simply related by reversing the direction of fermion flow in one of the two quark loops, and this does not change the expression for the diagram. The squared subprocess amplitudes, including an average over colour in the initial state, read Aqq¯→φ 2 = AΔqΔq¯→φ 2 = , Aδqδq¯→φ 2 = δjj′ c2Q2 , 2Nc where j, j′ are the indices for transverse quark polarisation and where we have used the spin projection operators given in equation (2.90) of [5]. Inserting (9.20) and (9.21) into (9.19), we obtain (9.21) (9.22) (9.23) dY dβ dσ1v1,pt = xx¯ g(x) g(x¯) 1 128π 1 − β2 Q4 (1 + β4) c4αs2 Z ∞ dy2 Nc2 − 1 0 y4 Φ2(yν) . Note that both (9.16) and (9.22) contain a product of gluon PDFs evaluated at the same x values, g(x) g(x¯). For the comparison, we can divide the common PDF factor out of the two expressions, in order to avoid having to use an explicit parameterisation. We also divide out various factors appearing in both expressions and compare the quantity Σ(β) = dσ Q2 dY dβ xx¯ g(x) g(x¯) 128π(Nc2 − 1) c4αs2 . It is straightforward to show that Σ is a function of the variable β only. We plot Σ for SPS and the subtraction term in figure 20, where the curve for Σ1v1,pt corresponds to the central choice ν = Q of the cutoff scale. We see that Σ1v1,pt(β) is indeed generally of the same order of magnitude as ΣSPS(β). The agreement is perhaps surprisingly good, given that Σ1v1,pt(β) involves integrating a low q approximation to the matrix element squared outside the region where the approximation is valid. The agreement gets worse towards the endpoints β → 0 and β → 1, which correspond to the high energy limit and the threshold limit, respectively. It is to be expected that the agreement is especially bad at these points, since in the subtraction term we effectively assume that the integration over the squared transverse momentum q2 of the scalar particles goes from zero to values of order Q2. For both β → 0 and β → 1 this assumption becomes a poor one: near threshold the phase space in q2 shrinks to zero, whilst in the high energy limit, q2 can go up to size sˆ ≫ Q2. 10 Summary Consistently incorporating the perturbative splitting of one parton into two is a highly nontrivial problem for the theoretical description of double parton scattering. DPS graphs SPS (ΣSPS), and the subtraction term (Σ1v1,pt). in which such splittings occur in both protons (1v1 graphs) overlap with loop corrections to single parton scattering. Another type of graph, typically referred to as 2v1, in which one parton pair arises from a perturbative splitting, and the other pair is an “intrinsic” one already existing at the nonperturbative level, overlaps with twist-four contributions to the cross section. Finally there is an overlap between DPS contributions where a splitting occurs in both protons only in the amplitude or its conjugate, and SPS/DPS interference graphs. We have presented a scheme to compute DPS and to consistently merge its contribution to the cross section with SPS and the other terms just mentioned. The scheme works in a similar manner for collinear and for TMD factorisation. Ultraviolet divergences that arise from perturbative splitting in a naive treatment of DPS are regulated by a cutoff function Φ(yν) in transverse position space. This avoids modification of the position space DPDs, which are defined via operator matrix elements in close analogy to single-parton distributions. In collinear factorisation, these DPDs hence evolve according to a homogeneous DGLAP equation, whilst their TMD counterparts satisfy a generalisation of the renormalisation group equations for single-parton TMDs. No modification of hard scattering cross sections computed for standard collinear or TMD factorisation is necessary in our scheme. Collins-Soper type equations describe the rapidity evolution of transversemomentum dependent DPDs and of collinear DPDs in colour non-singlet channels [21]. The problem of double counting between DPS and other contributions — notably between DPS and SPS — is solved by subtraction terms as specified in (4.16) and (4.22), which are obtained in a simple way from σDPS by replacing the DPDs with their appropriate short-distance limits. This paves the way for using the scheme at higher orders in αs, with calculations being considerably simpler for the subtraction terms compared with the full hard scattering process at the corresponding order. With a suitable choice of starting conditions and scales, specified in section 6.1, the DPS part of the cross section correctly resums DGLAP logarithms that are not included in the fixed order – 70 – Our scheme is naturally formulated with position-space DPDs F (xi, y), but it is possible to relate the Fourier transform of F (xi, y) Φ(yν) to DPDs F (xi, Δ) that are renormalised in transverse momentum space and satisfy an inhomogeneous DGLAP equation rather than a homogeneous one. This relation has the form of a perturbative matching equation, see (7.15), and is somewhat similar to the matching between PDFs defined in different schemes such as the MS and the DIS scheme. The momentum space representation also allows us to show that for the 2v1 contribution to DPS our scheme is equivalent to the ones in [12, 13] and in [14, 15] to leading logarithmic accuracy. For collinear DPDs, one can make a model ansatz consisting of two terms which in the limit y ≪ 1/Λ respectively give the perturbative splitting and the intrinsic part of the distribution. With such an ansatz, the DPS cross section naturally splits into 1v1, 2v1 and 2v2 terms, where 2v2 refers to contributions in which the parton pairs from both protons are intrinsic. A crucial question is how large DPS is compared with SPS at the perturbative order where graphs contribute to both mechanisms. This is especially acute in collinear factorisation, where DPS is power suppressed with respect to SPS. Note that only in very few channels (notably pair production of electroweak gauge bosons) SPS calculations are available at the required order. We argue in section 6.2 that in our scheme the variation with ν of the 1v1 term in DPS provides an order-of magnitude estimate for the SPS contribution σSPS (at the appropriate perturbative order), as it involves the same PDFs, overall coupling constants and kinematic region (small y, corresponding to large transverse momenta and virtualities of internal lines). An alternative estimate is provided by the double counting subtraction term σ1v1,pt, which by construction is dominated by small y. For the hypothetical process of scalar boson pair production from two gluons, we have shown that the latter estimate works well within about a factor of two, provided that one stays away from the extreme kinematic regions where the velocity β of a boson in the boson pair rest frame is close to 0 or 1. We constructed explicit (collinear) DPD input forms using the model ansatz just described, restricting ourselves to three quark flavours for simplicity and ease of implementation. These inputs were then numerically evolved to other scales using a code that implements the homogeneous double DGLAP equation. We used the resultant DPDs to compute so-called DPS luminosities (DPS cross sections omitting the process-dependent hard parts) and plotted these under various conditions. We observed that generically, the 1v1 contributions to the luminosity (both unpolarised and polarised) are comparable to or larger than the 2v1 and 2v2 contributions, with a large dependence on the cutoff ν. This demonstrates that, when including DPS in a cross section calculation, one must in general include σSPS up to the order that contains DPS-like double box graphs, together with the associated subtraction term (with unpolarised and polarised partons). Otherwise, one would have an uncertainty on the overall cross section that is as large as, or larger than the DPS term itself. We also confirmed that the ν variation of the 1v1 DPS contribution is indeed comparable to the central value of the associated double counting subtraction term, so that either of them may be used as an estimate for the SPS contribution. We identified several processes and scenarios where the ν variation of the 1v1 DPS luminosity is considerably smaller than the central value. As we argued above, one may HJEP06(217)83 contribution is reduced are when √ then justifiably neglect the order of σSPS containing the first DPS-like double box, as well as the associated subtraction term, compared with σDPS. One scenario of this kind is when the flavour indices in both DPDs are ud¯ (the luminosity with this parton flavour combination appears in W +W + production). The suppression of the DPS-like double box in σSPS is in this case related to the fact that perturbative splitting in the ud¯ DPD starts at order αs2 rather than at order αs. Further scenarios in which the ν variation in the 1v1 s becomes very large compared to the hard scales Qi, or when the rapidity separation between the produced hard systems is large. Both of these scenarios involve small x values in the DPDs — in the first, both x values in each DPD are small, whilst in the second, one x value in each DPD becomes much smaller than the other. Such processes and kinematic regions are the most promising ones to make useful calculations and measurements for DPS. In fact, several measurements investigating DPS have already been made in kinematics with Qi ≪ s or with large rapidity differences. It will be interesting in future work to make more complete and comprehensive predictions for such processes and kinematic regions in our framework, including for instance the full flavour dependence and contributions from all partonic channels for a considered final state. Fourier integrals sions. In this appendix we collect a number of results for Fourier transforms in 2 or 2 − 2ǫ dimen Fourier transform of a fractional power. To compute the Fourier transform of (k2)−λ in 2 − 2ǫ dimensions we write (k2)−λ = Γ−1(λ) R0∞ dα αλ−1e−αk2. One can then perform the Gaussian integration over k, and subsequently the integral over α. The result is Z d2−2ǫk e−iyk (k2)−λ = Γ(1 − ǫ − λ) π1−ǫ y2 −1+ǫ+λ As a corollary we obtain using kj eiyk = −i∂/(∂yj) eiyk. Z d2−2ǫk eiyk k j k2 = i 2 Γ(1 − ǫ) (4π)1−ǫ yj y2−2ǫ Integrals involving Φ. We now compute the integrals in (4.19), which are also needed for (7.6) and (7.7), If Φ is the step function Θ(u − b0) then 2 b0 J0(ur) = log r2 + b02 r2 1 4 2F3 1, 1; 2, 2, 2; − 4 b20 r2 where b0 is defined in (4.3) and 2F3(1, 1; 2, 2, 2; z) is a generalised hypergeometric function, which can be expanded as 1 + O(z) for small arguments. For Φ(u) = 1 − e−u2/4 we can proceed as follows: 2 0 u J0(ur)h1 − e−u2/4i = 2 = E1(r2) , Z 1/4 0 J0(ur) u2 dτ e−τu2 = Z 1/4 dτ 0 τ exp − 4τ (A.1) (A.2) (A.3) r 2 (A.4) γ + O(r2). The corresponding integrals with J2 instead of J0 are given by where E1 is the exponential integral. For small arguments one has E1(r2) = log(1/r2) − 2 b0 J2(ur) = 2J1(b0 r) b0 r and where in the second case we have used integration by parts. Both (A.5) and (A.6) behave like 1 + O(r2) for small r. Connection between the Fourier-Bessel transform and a cutoff. Consider the integral (8.7) with one term of the series (8.9). We will show that for Δ ≪ ν Z ∞ dy b0/ν y J0(yΔ) logn yν Z b0/Δ dy + O y Δ2 logn yν b0 + O logn−2 ν Δ . The following argument is similar to the derivation given in appendix B of [55]. The integral on the r.h.s. is readily performed and gives 2 J2(ur)h1 − e−u2/4i = Z ∞ 0 du J1(ur) ur ue−u2/4 = 1 − e−r2 , (A.6) (A.8) (A.10) Z ν/Δ 1 dv logn v = logn+1 ν . Δ The expression on the l.h.s. of (A.7) can be rewritten as Z ∞ b0Δ/ν du J0(u) logn uν b0Δ Z ∞ n + 1 b0Δ/ν du J1(u) logn+1 uν b0Δ +O Δ2 ν2 using integration by parts. Since J1(u) = O(u) at small u, the integral on the r.h.s. can be extended down to u = 0 with an accuracy of Δ2/ν2. Rewriting the logarithm using the binomial series, we obtain 1 X n+1 n + 1 n + 1 k=0 k logn+1 ν Δ logn+1−k ν Δ + X Z ∞ 0 k n+1 n + 1 du J1(u) logk u b0 dk logn+1−k ν Δ . In the last step we performed the integrals for k = 0, 1, 2 explicitly (the ones for k = 1, 2 are zero) and replaced the remaining ones by coefficients dk whose values are not important here. Comparison with (A.8) gives the desired result (A.7). Acknowledgments We gratefully acknowledge discussions with V. Braun, K. Golec-Biernat, D. Boer, A. Sch¨afer and F. Tackmann. Special thanks are due to T. Kasemets for valuable remarks on the manuscript. J.G. acknowledges financial support from the European Community under the FP7 Ideas program QWORK (contract 320389). Two of us (M.D. and J.G.) thank the Erwin Schr¨odinger International Institute for Mathematics and Physics (ESI) for hospitality during the programme “Challenges and Concepts for Field Theory and Applications in the Era of LHC Run-2”, when portions of this work were completed. The Feynman diagrams in this paper were produced with JaxoDraw [56, 57]. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. [1] R. Astalos et al., Proceedings of the Sixth International Workshop on Multiple Partonic Interactions at the Large Hadron Collider, arXiv:1506.05829. [2] H. Jung, D. Treleani, M. Strikman and N. van Buuren eds. Proceedings of the 7th International Workshop on Multiple Partonic Interactions at the LHC, DESY-PROC-2016-01, [INSPIRE]. [3] N. Paver and D. Treleani, Multi-Quark Scattering and Large pT Jet Production in Hadronic Collisions, Nuovo Cim. A 70 (1982) 215 [INSPIRE]. [4] M. Diehl and A. Sch¨afer, Theoretical considerations on multiparton interactions in QCD, Phys. Lett. B 698 (2011) 389 [arXiv:1102.3081] [INSPIRE]. [5] M. Diehl, D. Ostermeier and A. Sch¨afer, Elements of a theory for multiparton interactions in QCD, JHEP 03 (2012) 089 [Erratum ibid. 1603 (2016) 001] [arXiv:1111.0910] [INSPIRE]. [6] R. Kirschner, Generalized Lipatov-Altarelli-Parisi Equations and Jet Calculus Rules, [7] V.P. Shelest, A.M. Snigirev and G.M. Zinovev, The Multiparton Distribution Equations in Phys. Lett. B 84 (1979) 266 [INSPIRE]. QCD, Phys. Lett. B 113 (1982) 325 [INSPIRE]. [8] A.M. Snigirev, Double parton distributions in the leading logarithm approximation of perturbative QCD, Phys. Rev. D 68 (2003) 114012 [hep-ph/0304172] [INSPIRE]. [9] J.R. Gaunt and W.J. Stirling, Double Parton Distributions Incorporating Perturbative QCD Evolution and Momentum and Quark Number Sum Rules, JHEP 03 (2010) 005 [arXiv:0910.4347] [INSPIRE]. [10] F.A. Ceccopieri, An update on the evolution of double parton distributions, Phys. Lett. B 697 (2011) 482 [arXiv:1011.6586] [INSPIRE]. [11] M. Cacciari, G.P. Salam and S. Sapeta, On the characterisation of the underlying event, JHEP 04 (2010) 065 [arXiv:0912.4926] [INSPIRE]. [12] B. Blok, Yu. Dokshitser, L. Frankfurt and M. Strikman, pQCD physics of multiparton interactions, Eur. Phys. J. C 72 (2012) 1963 [arXiv:1106.5533] [INSPIRE]. [14] M.G. Ryskin and A.M. Snigirev, A fresh look at double parton scattering, [15] M.G. Ryskin and A.M. Snigirev, Double parton scattering in double logarithm approximation [16] A.V. Manohar and W.J. Waalewijn, What is Double Parton Scattering?, Phys. Lett. B 713 (2012) 196 [arXiv:1202.5034] [INSPIRE]. [17] J.R. Gaunt, Single Perturbative Splitting Diagrams in Double Parton Scattering, JHEP 01 (2013) 042 [arXiv:1207.0480] [INSPIRE]. [18] J. Collins, Foundations of perturbative QCD, Cambridge University Press, (2013). [19] T.C. Rogers, An Overview of Transverse Momentum Dependent Factorization and Evolution, Eur. Phys. J. A 52 (2016) 153 [arXiv:1509.04766] [INSPIRE]. [20] M. Diehl, J.R. Gaunt, D. Ostermeier, P. Plo¨ßl and A. Sch¨afer, Cancellation of Glauber gluon exchange in the double Drell-Yan process, JHEP 01 (2016) 076 [arXiv:1510.08696] [INSPIRE]. [21] M. Buffing, M. Diehl and T. Kasemets, Transverse momentum in double parton scattering: factorisation, evolution and matching, preprint in preparation DESY-2017-014, NIKHEF-2016-028 (2017). [22] A.V. Manohar and W.J. Waalewijn, A QCD Analysis of Double Parton Scattering: Color Correlations, Interference Effects and Evolution, Phys. Rev. D 85 (2012) 114009 [arXiv:1202.3794] [INSPIRE]. [23] H.-M. Chang, A.V. Manohar and W.J. Waalewijn, Double Parton Correlations in the Bag Model, Phys. Rev. D 87 (2013) 034009 [arXiv:1211.3132] [INSPIRE]. [24] M. Rinaldi, S. Scopetta and V. Vento, Double parton correlations in constituent quark models, Phys. Rev. D 87 (2013) 114021 [arXiv:1302.6462] [INSPIRE]. [25] M. Rinaldi, S. Scopetta, M. Traini and V. Vento, Double parton correlations and constituent quark models: a Light Front approach to the valence sector, JHEP 12 (2014) 028 [arXiv:1409.1500] [INSPIRE]. [26] M. Rinaldi, S. Scopetta, M.C. Traini and V. Vento, Correlations in Double Parton Distributions: Perturbative and Non-Perturbative effects, JHEP 10 (2016) 063 [arXiv:1608.02521] [INSPIRE]. [27] M. Rinaldi and F.A. Ceccopieri, Relativistic effects in model calculations of double parton distribution function, Phys. Rev. D 95 (2017) 034040 [arXiv:1611.04793] [INSPIRE]. [28] W. Broniowski and E. Ruiz Arriola, Valence double parton distributions of the nucleon in a simple model, Few Body Syst. 55 (2014) 381 [arXiv:1310.8419] [INSPIRE]. [29] W. Broniowski, E. Ruiz Arriola and K. Golec-Biernat, Generalized Valon Model for Double Parton Distributions, Few Body Syst. 57 (2016) 405 [arXiv:1602.00254] [INSPIRE]. [30] T. Kasemets and A. Mukherjee, quark-gluon double parton distributions in the light-front dressed quark model, Phys. Rev. D 94 (2016) 074029 [arXiv:1606.05686] [INSPIRE]. [31] V.M. Braun, A.N. Manashov and J. Rohrwild, Baryon Operators of Higher Twist in QCD and Nucleon Distribution Amplitudes, Nucl. Phys. B 807 (2009) 89 [arXiv:0806.2531] [32] V.M. Braun, A.N. Manashov and J. Rohrwild, Renormalization of Twist-Four Operators in QCD, Nucl. Phys. B 826 (2010) 235 [arXiv:0908.1684] [INSPIRE]. Nucl. Phys. B 894 (2015) 161 [arXiv:1405.2828] [INSPIRE]. [34] A.P. Bukhvostov, G.V. Frolov, L.N. Lipatov and E.A. Kuraev, Evolution Equations for Quasi-Partonic Operators, Nucl. Phys. B 258 (1985) 601 [INSPIRE]. [35] R.L. Jaffe, Spin, twist and hadron structure in deep inelastic processes, in The spin structure of the nucleon. Proceedings, International School of Nucleon Structure, 1st Course, Erice, Italy, August 3-10, 1995, pp. 42–129, (1996), hep-ph/9602236 [INSPIRE]. [36] P.J. Mulders and J. Rodrigues, Transverse momentum dependence in gluon distribution and fragmentation functions, Phys. Rev. D 63 (2001) 094021 [hep-ph/0009343] [INSPIRE]. [37] R.L. Jaffe and X.-D. Ji, Chiral odd parton distributions and Drell-Yan processes, Nucl. Phys. B 375 (1992) 527 [INSPIRE]. [38] J.C. Collins, D.E. Soper and G.F. Sterman, Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production, Nucl. Phys. B 250 (1985) 199 [INSPIRE]. [39] M. Diehl, T. Kasemets and S. Keane, Correlations in double parton distributions: effects of evolution, JHEP 05 (2014) 118 [arXiv:1401.1233] [INSPIRE]. [40] F. Caola, M. Dowling, K. Melnikov, R. R¨ontsch and L. Tancredi, QCD corrections to vector boson pair production in gluon fusion including interference effects with off-shell Higgs at the LHC, JHEP 07 (2016) 087 [arXiv:1605.04610] [INSPIRE]. [41] A.V. Belitsky and A.V. Radyushkin, Unraveling hadron structure with generalized parton distributions, Phys. Rept. 418 (2005) 1 [hep-ph/0504030] [INSPIRE]. [42] M. Diehl and J.R. Gaunt, Double parton scattering in the ultraviolet: addressing the double counting problem, in Proceedings, 7th International Workshop on Multiple Partonic Interactions at the LHC (MPI@LHC 2015): Miramare, Trieste, Italy, November 23–27, 2015, pp. 121–125, (2016), arXiv:1603.05468 [INSPIRE]. [43] M. Diehl and J.R. Gaunt, Double parton scattering in the ultraviolet: addressing the double counting problem, PoS(QCDEV2016)014 [arXiv:1611.01316] [INSPIRE]. [44] Z. Ligeti, I.W. Stewart and F.J. Tackmann, Treating the b quark distribution function with reliable uncertainties, Phys. Rev. D 78 (2008) 114014 [arXiv:0807.1926] [INSPIRE]. [45] R. Abbate, M. Fickinger, A.H. Hoang, V. Mateu and I.W. Stewart, Thrust at N3LL with Power Corrections and a Precision Global Fit for alphas(mZ), Phys. Rev. D 83 (2011) 074021 [arXiv:1006.3080] [INSPIRE]. [46] I.W. Stewart, F.J. Tackmann, J.R. Walsh and S. Zuberi, Jet pT resummation in Higgs production at N N LL′ + N N LO, Phys. Rev. D 89 (2014) 054001 [arXiv:1307.1808] [47] M. Mekhfi and X. Artru, Sudakov Suppression of Color Correlations in Multiparton Scattering, Phys. Rev. D 37 (1988) 2618 [INSPIRE]. K-Transverse, Nucl. Phys. B 197 (1982) 446 [INSPIRE]. transverse-momentum-dependent parton densities and the Collins-Soper evolution kernel, Phys. Rev. D 91 (2015) 074020 [arXiv:1412.3820] [INSPIRE]. contributions to double parton scattering production of two quarkonia, two Higgs bosons and cc¯cc¯, Phys. Rev. D 90 (2014) 054017 [arXiv:1407.5821] [INSPIRE]. JHEP 05 (2013) 150 [arXiv:1303.0842] [INSPIRE]. JHEP 01 (2013) 121 [arXiv:1210.5434] [INSPIRE]. Nucl. Phys. B 309 (1988) 282 [INSPIRE]. descriptions of semi-inclusive processes at low and high transverse momentum, JHEP 08 (2008) 023 [arXiv:0803.0227] [INSPIRE]. drawing Feynman diagrams. Version 2.0 release notes, Comput. Phys. Commun. 180 (2009) 1709 [arXiv:0811.4113] [INSPIRE]. multi-parton collisions , Eur. Phys. J. C 74 ( 2014 ) 2926 [arXiv: 1306 .3763] [INSPIRE]. of perturbative QCD , Phys. Rev. D 86 ( 2012 ) 014018 [arXiv: 1203 .2330] [INSPIRE]. [48] A.D. Martin , W.J. Stirling , R.S. Thorne and G. Watt , Parton distributions for the LHC , Eur. Phys. J. C 63 ( 2009 ) 189 [arXiv: 0901 .0002] [INSPIRE]. [49] J.C. Collins and D.E. Soper , Back- To-Back Jets: Fourier Transform from B to [50] J. Collins and T. Rogers , Understanding the large-distance behavior of [51] J.R. Gaunt , R. Maciula and A. Szczurek , Conventional versus single-ladder-splitting

This is a preview of a remote PDF:

Markus Diehl, Jonathan R. Gaunt, Kay Schönwald. Double hard scattering without double counting, Journal of High Energy Physics, 2017, 83, DOI: 10.1007/JHEP06(2017)083