Double hard scattering without double counting

TL;DR

This paper addresses the double counting issue in double parton scattering in proton-proton collisions by proposing a new operator-based framework, enabling higher-order calculations and better accuracy in distinguishing between single and double scattering effects.

Contribution

The authors introduce a novel operator-based method to define double parton distributions, resolving double counting issues and facilitating higher-order perturbative calculations in double parton scattering.

Findings

The framework allows estimation of higher-order contributions to single scattering.

Identification of conditions where higher-order effects are significant.

Guidance on when higher-order calculations are necessary for accuracy.

Abstract

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.

Tables

Table 1. Table 1: Examples for the cutoff function Φ Φ \Phi and relevant integrals. Exact results for the integral in the second row are given in ( 167 ) and ( A ). The integral in the third row converges for ϵ < 0 italic-ϵ 0 \epsilon<0 .

Equations

$$F_{\text{bare}}(x_i,{\mathbf{z}}_i,\mathbf{y})$$

$${\mathcal{O}}_i(y,z_i)$$

$$\frac{d}{d\log {\mu }_1^2}{F}_{a_1{a}_2}(x_1,x_2,\mathbf{y};{\mu }_1,{\mu }_2)$$

$$\displaystyle\qquad=\sum_{b_{1}}\int_{x_{1}}^{1-x_{2}}\frac{dx_{1}^{\prime}}{x_{1}^{\prime}}\;P_{a_{1}b_{1}}\Bigl{(}\frac{x_{1}}{x_{1}^{\prime}},\alpha_{s}(\mu_{1})\Bigr{)}\,F_{b_{1}a_{2}}(x_{1}^{\prime},x_{2},\boldsymbol{y};\mu_{1},\mu_{2})$$

$$x_i$$

$$d^4{q}_i$$

$$\frac{d{\sigma }_{\text{DPS}}}{d{x}_{1}d{x}_{2}d{\bar{x}}_{1}d{\bar{x}}_2{d}^2{\mathbf{q}}_1{d}^2{\mathbf{q}}_2}=\frac{1}{C}\sum _{a_1{a}_2{b}_1{b}_2}{\hat{\sigma }}_{a_1{b}_1}(Q_1^2,{\mu }_1^2){\hat{\sigma }}_{a_2{b}_2}(Q_2^2,{\mu }_2^2)$$

$$\times \int \frac{d^2{\mathbf{z}}_1}{(2\pi {)}^2}\frac{d^2{\mathbf{z}}_2}{(2\pi {)}^2}{d}^2\mathbf{y}{e}^{-i({\mathbf{q}}_1{\mathbf{z}}_1+{\mathbf{q}}_2{\mathbf{z}}_2)}{F}_{b_1{b}_2}({\bar{x}}_i,{\mathbf{z}}_i,\mathbf{y};{\mu }_i)F_{a_1{a}_2}(x_i,{\mathbf{z}}_i,\mathbf{y};{\mu }_i)$$

$$\frac{d{\sigma }_{\text{DPS}}}{d{x}_{1}d{x}_{2}d{\bar{x}}_{1}d{\bar{x}}_2}=\frac{1}{C}\sum _{a_1{a}_2{b}_1{b}_2}{\int }_{x_1}^{1-x_2}\frac{d{x}_1^{\prime }}{x_1^{\prime }}{\int }_{x_2}^{1-x_1^{\prime }}\frac{d{x}_2^{\prime }}{x_2^{\prime }}{\int }_{{\bar{x}}_1}^{1-{\bar{x}}_2}\frac{d{\bar{x}}_1^{\prime }}{{\bar{x}}_1^{\prime }}{\int }_{{\bar{x}}_2}^{1-{\bar{x}}_1^{\prime }}\frac{d{\bar{x}}_2^{\prime }}{{\bar{x}}_2^{\prime }}$$

$$\times {\hat{\sigma }}_{a_1{b}_1}(x_1^{\prime }{\bar{x}}_1^{\prime }s,{\mu }_1^2){\hat{\sigma }}_{a_2{b}_2}(x_2^{\prime }{\bar{x}}_2^{\prime }s,{\mu }_2^2)\int d^2\mathbf{y}{F}_{b_1{b}_2}({\bar{x}}_i^{\prime },\mathbf{y};{\mu }_i)F_{a_1{a}_2}(x_i^{\prime },\mathbf{y};{\mu }_i)$$

$$\frac{1}{(2\pi {)}^{4-4ϵ}}\int d^{2-2ϵ}{\mathbf{z}}_1{d}^{2-2ϵ}{\mathbf{z}}_2{e}^{-i({\mathbf{z}}_1{\mathbf{k}}_1+{\mathbf{z}}_2{\mathbf{k}}_2)}$$

$$\int d^{2-2ϵ}\mathbf{y}{e}^{i\mathbf{y}\mathbf{\Delta }}$$

$$\mathbf{q}$$

$$\frac{d\sigma }{d^2{\mathbf{q}}_1{d}^2{\mathbf{q}}_2}\sim \frac{1}{Q_i^4{q}_i^2}.$$

$${\mathbf{y}}_{\pm }$$

$${\mathbf{k}}_{\pm }$$

$$F(x_i,{\mathbf{y}}_{\pm },\mathbf{Z})$$

$$=\frac{1}{(2\pi {)}^{2-2ϵ}}\int d^{2-2ϵ}\mathbf{K}{d}^{2-2ϵ}{\mathbf{k}}_{+}{d}^{2-2ϵ}{\mathbf{k}}_{-}{e}^{i\mathbf{Z}\mathbf{K}+i({\mathbf{y}}_{+}{\mathbf{k}}_{-}-{\mathbf{y}}_{-}{\mathbf{k}}_{+})}F(x_i,{\mathbf{k}}_{\pm },\mathbf{K})$$

$$F_{a_1{a}_2,\text{spl,pt}}(x_i,{\mathbf{k}}_{\pm },\mathbf{K})$$

$$\displaystyle\quad\times\biggl{[}\frac{f_{a_{0}}(x_{1}+x_{2},\boldsymbol{K})}{x_{1}+x_{2}}\,\frac{\alpha_{s}}{2\pi}\,T^{jj^{\prime}}_{a_{0}\to a_{1}a_{2}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}},\epsilon\biggr{)}+\ldots\,\biggr{]}\,,$$

$$F_{a_1{a}_2,\text{spl,pt}}(x_i,{\mathbf{y}}_{\pm },\mathbf{Z})$$

$$\displaystyle\quad\times\biggl{[}\frac{f_{a_{0}}(x_{1}+x_{2},\boldsymbol{Z})}{x_{1}+x_{2}}\,\frac{\alpha_{s}}{2\pi}\,T^{jj^{\prime}}_{a_{0}\to a_{1}a_{2}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}},\epsilon\biggr{)}+\ldots\,\biggr{]}$$

$$\int d^2\mathbf{y}{d}^2{\mathbf{z}}_1{d}^2{\mathbf{z}}_2{e}^{-i({\mathbf{q}}_1{\mathbf{z}}_1+{\mathbf{q}}_2{\mathbf{z}}_2)}$$

$$\frac{1}{\pi }\int d^2\mathbf{y}{e}^{\pm i\mathbf{q}\mathbf{y}}\frac{{\mathbf{y}}^j{\mathbf{y}}^{j^{\prime }}}{y^4}$$

$$F_{{\alpha }_1{\alpha }_2,y_{-}\to 0}(x_i,{\mathbf{y}}_{\pm },\mathbf{Z})$$

$$\displaystyle\bigl{\langle}p\mskip 1.5mu\big{|}\phi_{0}\bigl{(}-{\tfrac{1}{2}}Z\bigr{)}\,\phi_{2}\bigl{(}{\tfrac{1}{2}}(Z-y_{+})\bigr{)}\,\phi_{1}\bigl{(}{\tfrac{1}{2}}(Z+y_{+})\bigr{)}\big{|}\mskip 1.5mup\bigr{\rangle}$$

$$\displaystyle\qquad\quad=\bigl{\langle}p\mskip 1.5mu\big{|}\phi_{0}\bigl{(}{\tfrac{1}{2}}(y-Z)\bigr{)}\,\phi_{2}\bigl{(}{\tfrac{1}{2}}z_{2}\bigr{)}\,\phi_{1}\bigl{(}y+{\tfrac{1}{2}}z_{1}\bigr{)}\big{|}\mskip 1.5mup\bigr{\rangle}\,,$$

$$D_{{\alpha }_1{\alpha }_2|{\alpha }_0,y_{+}\to 0}(x_i,{\mathbf{y}}_{+},\mathbf{Z})$$

$$F_{{\alpha }_1{\alpha }_2,y_{\pm }\to 0}(x_i,{\mathbf{y}}_{\pm },\mathbf{Z})$$

$$\displaystyle\quad=\frac{\boldsymbol{y}_{+}^{j}}{y_{+}^{2}}\,\frac{\boldsymbol{y}_{-}^{j^{\prime}}}{y_{-}^{2}}\,U^{j}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}(x_{i})\,\bigl{[}U^{j^{\prime}}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}(x_{i})\bigr{]}^{*}\,f_{\alpha_{0}}(x_{1}+x_{2},\boldsymbol{Z})\,.$$

$$F_{a_1{a}_2,\text{spl,pt}}(x_1,x_2,\mathbf{y};\mu )$$

$$\displaystyle\quad=\frac{\mu^{2\epsilon}}{y^{2-4\epsilon}}\,\frac{\Gamma^{2}(1-\epsilon)}{\pi^{1-2\epsilon}}\,\frac{f_{a_{0}}(x_{1}+x_{2};\mu)}{x_{1}+x_{2}}\,\frac{\alpha_{s}(\mu)}{2\pi}\,P_{a_{0}\to a_{1}a_{2}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}},\epsilon\biggr{)}\,.$$

$$P_{g\to q\bar{q}}(u,ϵ)$$

$${}^1{F}_{a_1{a}_2,\text{spl,pt}}(x_1,x_2,\mathbf{y};\mu )$$

$$\displaystyle\qquad=\frac{1}{\pi^{1-\epsilon}\,y^{2-2\epsilon}}\,\sum_{a_{0}}\int\limits_{x_{1}+x_{2}}^{1}\frac{dv}{v}\,\frac{f_{a_{0}}(v;\mu)}{v}\,V_{a_{0}\to a_{1}a_{2}}\biggl{(}\frac{x_{1}}{v},\frac{x_{2}}{v},\alpha_{s}(\mu),y\mu,\epsilon\biggr{)}$$

$$V_{a_0\to a_1{a}_2}(v_1,v_2,{\alpha }_s,y\mu ,ϵ)$$

Full text

1 Introduction

The precise description of high-energy proton-proton collisions in QCD is imperative for maximising the physics potential of the LHC and of possible future hadron colliders. An important issue in this context is to understand the mechanism of double parton scattering (DPS), in which two pairs of partons undergo a hard scattering in one and the same proton-proton collision. In many situations DPS is suppressed compared with single parton scattering (SPS), but this suppression generically becomes weaker with increasing collision energy. For specific kinematics or specific final states, DPS can become comparable to or even larger than SPS. An overview of recent experimental and theoretical activities in this area can for instance be found in Astalos:2015ivw ; Proceedings:2016tff .

Consider a DPS process $pp\to Y_{1}+Y_{2}+X$, where $Y_{1}$ and $Y_{2}$ are observed particles or groups of particles produced in two separate hard scattering processes, whilst $X$ denotes all unobserved particles in the final state. A cross section formula has been put forward long ago in the framework of collinear factorisation, where the transverse momenta $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ of $Y_{1}$ and $Y_{2}$ are integrated over Paver:1982yp . A corresponding expression has been given in Diehl:2011tt ; Diehl:2011yj for transverse momentum dependent (TMD) factorisation, where $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ are small compared with the hard scales in $Y_{1}$ and $Y_{2}$.

Inside a proton, the two partons that take part in the two hard scatters can originate from the perturbative splitting of one parton. The relevance of this splitting mechanism for the evolution equations of double parton distributions (DPDs) has been realised long ago Kirschner:1979im ; Shelest:1982dg and studied more recently in Snigirev:2003cq ; Gaunt:2009re ; Ceccopieri:2010kg . However, it was only noted in Diehl:2011tt ; Diehl:2011yj 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 Cacciari:2009dp .

Different ways of dealing with this issue have been proposed in Blok:2011bu ; Blok:2013bpa , Ryskin:2011kk ; Ryskin:2012qx and Manohar:2012pe . 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 $\alpha_{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 $\alpha_{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 $\boldsymbol{y}$ between the two partons, but we show in section 7 how one may instead use DPDs depending on the transverse momentum conjugate to $\boldsymbol{y}$. This allows us to compare our results with those of Blok:2011bu ; Blok:2013bpa and Ryskin:2011kk ; Ryskin:2012qx , 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.111We follow here the nomenclature of Gaunt:2012dd . The 2v1 contribution is referred to as $4\times 2$ in Diehl:2011yj , as $3\to 4$ in Blok:2011bu ; Blok:2013bpa (where four-jet production is considered), and as $1\times 2$ in Ryskin:2011kk ; Ryskin:2012qx . We also give some basic definitions and results. Throughout this work we use light-cone coordinates $v^{\pm}=(v^{0}+v^{3})/\sqrt{2}$ for any four-vector $v$. We write $\boldsymbol{v}=(v^{1},v^{2})$ for the transverse components and $v=\sqrt{\boldsymbol{v}^{2}}$ for the length of the transverse vector.

Since the perturbative splitting mechanism in DPDs leads to issues in the ultraviolet (UV) region, renormalisation plays a crucial role in our context. We work in dimensional regularisation and extend the definitions of Diehl:2011yj to $D=4-2\epsilon$ dimensions, changing phase space factors $(2\pi)^{4}$ into $(2\pi)^{D}$. Bare two-parton distributions are then given by

[TABLE]

with twist-two operators

[TABLE]

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 $\Gamma_{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 (1).

In the process of deriving factorisation, one finds that the proton matrix element in (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 Collins:2011zzd and in the recent overview Rogers:2015sqa . A brief account for DPDs is given in section 2.1 of Diehl:2015bca , and a more detailed discussion will be provided in Buffing:2016wip . The product $F_{\text{bare}}\times\{\text{soft factors}\}$ depends on a parameter $\zeta$ which regulates rapidity divergences. A scheme in which a soft factor appears explicitly in the cross section formula was presented in Manohar:2012jr . Since the treatment of soft gluons does not affect parton splitting in any special way, we will not discuss it further in the present paper. Correspondingly, we will suppress the argument $\zeta$ in all DPDs.

As a final step one performs UV renormalisation, which we assume to be done in the $\overline{\text{MS}}$ scheme. The DPDs obtained from (1) are appropriate for TMD factorisation, with the transverse positions $\boldsymbol{z}_{i}$ being Fourier conjugate to the transverse parton momenta that determine the final state kinematics. These distributions renormalise multiplicatively, with one renormalisation factor for each product of a quark field (or gluon field strength) with the Wilson line at the same transverse position and one factor for each pair of Wilson lines at the same transverse position in the soft factors. Denoting the product of these factors with $Z$, one obtains the final DPD as F=\lim_{\epsilon\to 0}\bigl{(}Z\times F_{\text{bare}}\times\{\text{soft factors}\}\bigr{)}. It obeys a renormalisation group equation which is a straightforward generalisation of the one for TMDs (given e.g. in  Rogers:2015sqa ).

The DPDs needed for collinear factorisation are obtained by setting $\boldsymbol{z}_{i}=\boldsymbol{0}$ in $F_{\text{bare}}$ and in the associated soft factors, before renormalisation. Setting $\boldsymbol{z}_{i}=\boldsymbol{0}$ introduces ultraviolet divergences in the operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, and in the associated soft factors. The renormalised DPDs are then obtained as F=\lim_{\epsilon\to 0}\bigl{(}Z_{1}\otimes Z_{2}\otimes F_{\text{bare}}\times\{\text{soft factors}\}\bigr{)}, where $Z_{i}$ renormalises the operators associated with parton $i$ and where the convolution products are in the momentum fractions $x_{i}$. In the colour singlet channel, where both operators $\mathcal{O}_{i}$ in (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 $\mu_{1}$ and $\mu_{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 $\mu_{1}$ and additional jet production at a much lower scale $\mu_{2}$ (of course $\mu_{2}$ needs to remain in the perturbative region for our factorisation approach to be justified).

With different scales $\mu_{1},\mu_{2}$ in the collinear DPDs, we have a homogeneous evolution equation

[TABLE]

in $\mu_{1}$ and its analogue for $\mu_{2}$. For colour singlet DPDs, the kernels $P_{a_{1}b_{1}}$ on the r.h.s. are the usual DGLAP splitting functions for single parton densities. In colour non-singlet channels, both the DPDs and the splitting kernels have an additional dependence on the rapidity parameter $\zeta$ Buffing:2016wip . Following the notation of Diehl:2011yj , the labels $a_{i},b_{i}$ in (2) denote both the species and the polarisation of the partons. The relevant labels are $q,\Delta q,\delta q$ for unpolarised, longitudinally polarised and transversely polarised quarks, likewise $\bar{q},\Delta\bar{q},\delta\bar{q}$ for antiquarks, and $g,\Delta g,\delta 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.

To simplify our presentation, we will consider the production $pp\to V_{1}+V_{2}+X$ of two electroweak gauge bosons $V_{i}=\gamma^{*},Z,W$. Our results readily generalise to other processes for which DPS factorisation can be established; in the case of TMD factorisation this requires that the produced particles are colour singlets. We denote the four-momenta of the two bosons by $q_{i}$, their squared invariant masses by $Q_{i}^{2}=q_{i}^{\mu}\mskip 1.5muq_{i\mu}$ and their rapidities by $Y_{i}={\tfrac{1}{2}}\log(q_{i}^{+}/q_{i}^{-})$. We work in the proton-proton centre-of-mass frame, taking the proton with momentum $p$ ($\bar{p}$) to move in the positive (negative) $3$ direction. Furthermore we define

[TABLE]

with $s=(p+\bar{p})^{2}$. For the phase space of each gauge boson we have

[TABLE]

The “naive” cross section formulae (not taking into account the UV problems discussed below) read

[TABLE]

for TMD factorisation and

[TABLE]

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 Diehl:2011yj , 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) or (2), and we do not treat them explicitly in the present work for ease of notation.

The hard scattering cross sections $\hat{\sigma}_{i}$ in (2) are for the exclusive final state $V_{i}$ with transverse momentum $\boldsymbol{q}_{i}$, which must satisfy $q_{i}\ll Q_{i}$. By contrast, $\hat{\sigma}_{i}$ in (2) is integrated over $\boldsymbol{q}_{i}$ and inclusive for final states $V_{i}+X$. At leading order in $\alpha_{s}$ it contains $\delta$ functions that enforce $x_{i}^{\prime}=x_{i}$ and $\bar{x}_{i}^{\prime}=\bar{x}_{i}$. The subtractions for collinear and soft regions in $\hat{\sigma}_{i}$ are different in the two factorisation frameworks, but in both cases they lead to a dependence on the factorisation scale $\mu_{i}$ that cancels against the $\mu_{i}$ dependence of the DPDs, up to powers in $\alpha_{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.

As was pointed out in Diehl:2011yj , the framework discussed so far suffers from problems in the region of small transverse distances between the two partons in a DPD. The leading behaviour of the collinear distributions $F(x_{i},\boldsymbol{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) this gives a quadratically divergent integral at small $y$, which clearly signals an inappropriate treatment of the ultraviolet region. As we will review in section 3.1, the short-distance behaviour of the distributions $F(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$ is less singular but still leads to logarithmic divergences in the TMD factorisation formula (2).

Instead of using DPDs depending on transverse positions, one may Fourier transform them to transverse momentum space, integrating

[TABLE]

for TMDs and collinear distributions. For collinear distributions, this transformation must be made before subtracting UV divergences and setting $\epsilon\to 0$: with $F(x_{i},\boldsymbol{y})\sim y^{-2}$ in $D=4$ dimensions, the Fourier integral (9) would be logarithmically divergent at $y=0$. In $D=4-2\epsilon$ dimensions this singularity turns into a pole in $1/\epsilon$, which requires an additional subtraction as we will review in section 7. Rather than being associated with the individual operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, this subtraction is related to the singularity arising when the transverse distance $y$ between the two operators goes to zero. It leads to an inhomogeneous evolution equation for the DPD $F(x_{i},\boldsymbol{\Delta})$ in momentum space, which has been extensively discussed in the literature Kirschner:1979im ; Shelest:1982dg ; Snigirev:2003cq ; Gaunt:2009re ; Ceccopieri:2010kg for the case $\boldsymbol{\Delta}=\boldsymbol{0}$. Notice, however, that this extra $\mu$ dependence does not cancel in the cross section when (2) is rewritten in transverse momentum space. Moreover, the additional UV renormalisation of $F(x_{i},\boldsymbol{\Delta})$ does not remove all UV divergences at the cross section level. The singularity of $F(x_{i},\boldsymbol{y})$ at small $y$ translates into a behaviour $F(x_{i},\boldsymbol{\Delta})\sim\log(\mu^{2}/\Delta^{2})$ at large $\Delta$ (see section 7), which gives a quadratic divergence for the $\boldsymbol{\Delta}$ integration in the DPS cross section.

It is easy to identify the origin of the UV divergences just discussed. Both in the $\boldsymbol{y}$ and $\boldsymbol{\Delta}$ 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 $\Delta\ll Q_{i}$ or, equivalently, $y\gg 1/Q_{i}$ (see section 2.1.2 in Diehl:2011yj ). 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 $\boldsymbol{q}_{1}$ and $\boldsymbol{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

[TABLE]

For $q\ll Q_{i}$ the graph gives a leading contribution to $d\sigma/d{q}^{2}$ if the transverse quark momenta in the loops are all of order $Q_{i}$. This carries the quark lines far off shell, so that this contribution is naturally associated with SPS, with $g+g\to V_{1}V_{2}$ as the hard subprocess. A leading contribution is also obtained from the region where all transverse quark momenta are much smaller than $Q_{i}$. This region is naturally described as DPS, with two disconnected hard scattering processes $q\bar{q}\to V_{1}$ and $q\bar{q}\to V_{2}$ and double parton distributions with perturbative $g\to q\bar{q}$ 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 $\boldsymbol{q}$.

Let us now turn to the graph in figure 1b and consider the cross section integrated over $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$. In the region of large $q\sim Q_{i}$, the quark lines at the top of the graph have transverse momenta of order $q$ and are far off shell. The proper description of this region is in terms of a hard scattering $q\bar{q}+g\to V_{1}V_{2}$, convoluted with a collinear twist-two 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\ll Q_{i}$, the $g\to q\bar{q}$ splittings are near collinear and the approximations for DPS are appropriate. We call this the 2v1 contribution to DPS, recalling that there is a $q\bar{q}+g$ subprocess in the graph. Both small and large $q$ give leading contributions to the integrated cross section, and in a naive calculation adding up the twist-four term and the DPS term has again a double counting problem, as well as divergences in each contribution. The naive DPS cross section has a logarithmic divergence at small $y$, which is seen by inserting the $1/y^{2}$ splitting behaviour of only one DPD in (2). In turn, the hard scattering cross section in the twist-four term contains a collinear divergence in the form of an integral behaving like $dq^{2}/q^{2}$ at $q\to 0$, as we will show in section 4.1.2.

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.222The intrinsic part of a DPD may be studied using quark models Chang:2012nw ; Rinaldi:2013vpa ; Rinaldi:2014ddl ; Rinaldi:2016jvu ; Rinaldi:2016mlk ; Broniowski:2013xba ; Broniowski:2016trx ; Kasemets:2016nio , 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. This has been pursued independently by Blok et al. Blok:2011bu ; Blok:2013bpa and by Ryskin and Snigirev Ryskin:2011kk ; Ryskin:2012qx . 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 Ryskin:2011kk ; Ryskin:2012qx , where an ultraviolet cutoff must be imposed to regulate the quadratic divergence we mentioned earlier. By contrast, the authors of Blok:2011bu ; Blok:2013bpa 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 Collins:2011zzd 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 Diehl:2011yj 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\to V_{1}+V_{2}+X$ as a concrete example, but the discussion readily generalises to other cases.

The differential cross section $d\sigma/(d^{2}\boldsymbol{q}_{1}\mskip 1.5mud^{2}\boldsymbol{q}_{2})$ in the region $q_{i}\ll Q_{i}$ can be computed using TMD factorisation. Here and in the following we write $Q_{i}$ to denote the generic size of $Q_{1}$ and $Q_{2}$, and likewise for $q_{i}$. The transverse momenta $q_{i}$ may be of nonperturbative size $\Lambda$ or much larger. In the latter case, further simplifications are possible by expressing transverse momentum dependent distributions in terms of collinear ones Buffing:2016wip , but we shall not discuss this here. The leading power behaviour of the cross section is

[TABLE]

When $q_{i}$ goes to zero, it should be replaced by $\Lambda$. Three types of mechanisms contribute to the leading behaviour, namely DPS, SPS, and the interference between SPS and DPS. Corresponding graphs are shown in figures 1c, 2a and 2b, respectively.

As discussed in the previous section, certain graphs contribute both to DPS and to SPS, depending on the kinematics of their internal lines. The 1v1 graph in figure 1 also has leading regions in which one of the loops is hard and the other is collinear. These regions contribute to the SPS/DPS interference, as shown in figure 2c. The double counting problem thus concerns both SPS, DPS and their interference. Note that the SPS graph in figure 2a contributes to the SPS/DPS interference but not to DPS.

Both the amplitude and its conjugate in the 1v1 graph contains a loop integral that behaves like $d^{2}\boldsymbol{k}/\boldsymbol{k}^{2}$ in the region $\Lambda,q_{i}\ll k\ll Q_{i}$. When integrated over the full phase space, each loop thus builds up a so-called DPS logarithm, namely $\log(Q_{i}/q_{i})$ when q_{i}\raisebox{-4.0pt}{,\stackrel{{\scriptstyle\textstyle>}}{{\sim}},}\Lambda and $\log(Q_{i}/\Lambda)$ when q_{i}\raisebox{-4.0pt}{,\stackrel{{\scriptstyle\textstyle<}}{{\sim}},}\Lambda. Whether these logarithms reside in SPS, DPS or their interference depends on how exactly one handles the double counting problem. We come back to this issue in section 4.2.

Let us now turn to the cross section integrated over $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$, which can be described using collinear factorisation. The leading power behaviour of the cross section, $\sigma\sim 1/Q_{i}^{2}$, is given by the SPS mechanism alone. Several mechanisms contribute at the power suppressed level $\sigma\sim\Lambda^{2}/Q_{i}^{4}$. These are

DPS, which is suppressed because it can only populate the region $q_{i}\ll Q_{i}$ rather than the full phase space up to $q_{i}\sim Q_{i}$, 2. 2.

the interference between SPS and DPS, which is suppressed for the same reason, 3. 3.

hard scattering with a twist-four distribution for one proton and a twist-two distribution for the other. Example graphs are figure 3a, as well as figure 1b with the box removed. 4. 4.

hard scattering with twist-three distributions for both protons. An example graph is figure 3b.

The rationale for considering such contributions is that – whilst being power suppressed compared with SPS – they may be enhanced by higher parton luminosities at small momentum fractions $x$, or by coupling constants in the relevant hard scattering subprocesses. Let us emphasise that a complete calculation of the cross section at the level of $\Lambda^{2}/Q_{i}^{2}$ 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 Diehl:2011yj ). 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 Braun:2008ia ; Braun:2009vc ; Ji:2014eta . Twist-$n$ operators with $n$ parton fields were called “‘quasipartonic” in Bukhvostov:1985rn and involve only the “good” parton fields in the parlance of light-cone quantisation Jaffe:1996zw . 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 $\pm 1$ cannot add up to zero. For quark-antiquark-gluon operators it forces the quark and antiquark fields to have opposite chirality, i.e. one only has the operator combination $\bar{q}\sigma^{+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 Mulders:2000sh , whereas the corresponding quark-antiquark distributions are again chiral-odd Jaffe:1991ra . 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\sigma/dq^{2}\sim 1/q^{2}$ in the region $\Lambda\ll q\ll Q_{i}$, which gives a $\log(Q_{i}/\Lambda)$ 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 Diehl:2011yj we give expressions in $D=4-2\epsilon$ 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

[TABLE]

with Fourier conjugate momenta 333The momentum $\boldsymbol{\Delta}$ is called $\boldsymbol{r}$ in Diehl:2011tt ; Diehl:2011yj .

[TABLE]

The relation between DPDs in position and momentum space reads

[TABLE]

in $D=4-2\epsilon$ dimensions. As seen in figure 4, one can identify $\boldsymbol{y}_{+}$ ($\boldsymbol{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 $\boldsymbol{k}_{-}$ ($\boldsymbol{k}_{+}$). The splitting singularities of the DPDs thus occur at $y_{\pm}\to 0$ or $k_{\pm}\to\infty$.

The perturbative splitting contribution $F_{\text{spl,pt}}$ to transverse-momentum dependent DPDs in momentum space has been calculated at leading order in section 5.2.2 of Diehl:2011yj . Generalising these results to $D=4-2\epsilon$ dimensions, we have

[TABLE]

where $j$, $j^{\prime}$ are transverse Lorentz indices and $f_{a_{0}}(x_{1}+x_{2},\boldsymbol{K})$ is an unpolarised single-parton TMD.444Compared with section 5.2 of Diehl:2011yj , the kernel $T^{jj^{\prime}}$ used here has the opposite order of indices $jj^{\prime}$ and includes a colour factor, e.g. $T_{F}=1/2$ for the colour singlet distribution ${}^{1}F_{g\to q\bar{q}}$. The ellipsis denotes a term that involves a TMD for polarised partons in an unpolarised proton and depends on $\boldsymbol{K}$ but not on $\boldsymbol{k}_{\pm}$. In position space we then get

[TABLE]

using the Fourier integral (166), where the term denoted by an ellipsis depends on $\boldsymbol{Z}$ but not on $\boldsymbol{y}_{\pm}$. It is understood that for transverse quark or linear gluon polarisation, both $F_{a_{1}a_{2}}$ and the kernel $T$ carry additional transverse Lorentz indices. $f_{a_{0}}(x_{1}+x_{2},\boldsymbol{Z})$ is the Fourier transform of $f_{a_{0}}(x_{1}+x_{2},\boldsymbol{K})$. The form (3.1) gives the leading behaviour of the DPD for large $k_{\pm}\gg\Lambda$, and correspondingly (3.1) gives the leading behaviour for $y_{\pm}\ll 1/\Lambda$.

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

[TABLE]

with $\boldsymbol{q}$ defined in (10). Performing the angular integration in

[TABLE]

where $\boldsymbol{y}$ stands for $\boldsymbol{y}_{+}$ or $\boldsymbol{y}_{-}$, we see that the integral with $J_{0}$ is divergent at $y=0$. Given the range of validity of (3.1) one should impose $y\ll 1/\Lambda$ in (18), although the integrals are finite for $y\to\infty$ due to the oscillations of the Bessel functions.

The perturbative splitting contribution to DPDs at higher order in $\alpha_{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_{\pm}=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

[TABLE]

were $U_{\alpha_{0}\to\alpha_{1}\alpha_{2}}$ is a kernel for the splitting $\alpha_{0}\to\alpha_{1}\alpha_{2}$ in the amplitude (hence its complex conjugate appears in (19)). $D_{\alpha_{1}\alpha_{2}|\alpha_{0}}$ is the position space version of a transverse-momentum dependent twist-three distribution, constructed from the hadronic matrix element

[TABLE]

where $\phi_{i}$ is a “good” field for parton $\alpha_{i}$ (cf. section 2.1). Distributions $D_{\alpha_{0}|\alpha_{1}\alpha_{2}}(x_{i},\boldsymbol{y}_{-},\boldsymbol{Z})$ where $\alpha_{0}$ belongs to the amplitude and $\alpha_{1},\alpha_{2}$ to the complex conjugate amplitude are defined in analogy. In the second step of (3.1) we have used translation invariance and shifted the parton fields to the same position as in the corresponding DPD (see figure 4). The labels $\alpha_{i}$ denote the parton species; it is understood that in (19) and the following equations parton helicities are taken fixed on the l.h.s. and must be appropriately summed over on the r.h.s. Note the difference between this notation and the labels $a_{i}$, which denote parton species and polarisation (none, longitudinal, transverse or linear) and thus refer to a pair of parton legs. The notation with $a_{i}$ is hence not suitable for distributions with three parton fields.

If $y_{+}$ is small, then $D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})$ itself can be generated by perturbative splitting, as shown in figure 5c. We have

[TABLE]

Notice that a quark and antiquark produced by perturbative splitting have opposite helicities, so that the corresponding quark-antiquark operator $\phi_{2}\,\phi_{1}$ in $D_{q\bar{q}\mskip 1.5mu|g}$ must be chiral-even. Inserting (21) into (19) we obtain

[TABLE]

Taking appropriate linear combinations of parton helicities, we recover the form of $F_{a_{1}a_{2},\mskip 1.5mu\text{spl,pt}}$ in (3.1) at $\epsilon=0$.

3.2 Collinear DPDs: splitting contribution

We now turn to collinear DPDs, i.e. to the case where $\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{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 Diehl:2011yj . For $4-2\epsilon$ dimensions, one finds the general form

[TABLE]

The kernel for the splitting $g\to q\bar{q}$ reads for instance

[TABLE]

with a factor $f=1$ for the colour singlet and $f=-1/\sqrt{N^{2}-1}$ for the colour octet DPD. In terms of the kernel in (3.1) we have $T^{jj^{\prime}}_{g\to q\bar{q}}=\delta^{jj^{\prime}}_{\phantom{q}}P_{g\to q\bar{q}}$. We recognise in $P_{g\to q\bar{q}}(u,0)$ the usual DGLAP splitting kernel without the terms proportional to $\delta(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

[TABLE]

in $D=4-2\epsilon$ 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 $\epsilon=0$. The splitting kernel $V$ is a double series

[TABLE]

The $\mu$ (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 $\mu^{2\epsilon}\alpha_{s}(\mu)$; terms with $n>m$ in (26) 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)}(v_{1},v_{2},\epsilon)=\delta(1-v_{1}-v_{2})\,V^{(1)}(v_{1},\epsilon)$. Inserted into (3.2) this gives a form consistent with (3.2). Using that at order $\alpha_{s}^{n}$ the poles of highest order are $1/\epsilon^{n-1}$, we find

[TABLE]

in the physical limit $\epsilon=0$.

For colour nonsinglet DPDs one must regulate rapidity divergences, which complicates the preceding result. Taking e.g. Wilson lines along non-lightlike paths introduces additional vectors and changes the analysis of boost properties of the kernel. We will not pursue this issue here.

DPDs with transverse quark or linear gluon polarisation carry transverse Lorentz indices. Their perturbative splitting expressions thus have a tensor structure containing additional factors of $\boldsymbol{y}^{j}/y$ compared with the formulae above. At leading order one readily finds from (3.1) and the appropriate splitting kernels that the factor $1/y^{2-4\epsilon}$ in (3.2) is to be replaced with $\boldsymbol{y}^{j}\boldsymbol{y}^{j^{\prime}}/y^{4-4\epsilon}$ 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

[TABLE]

omitting all arguments other than $\boldsymbol{y}$ and $\boldsymbol{\Delta}$. 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 $\boldsymbol{\Delta}$, so that one can replace $F(\boldsymbol{\Delta})$ by its approximation for large $\Delta$, following the power counting analysis of section 5.2 in Diehl:2011yj . This leads us to write

[TABLE]

where the three terms on the r.h.s. will be described shortly. In $D=4$ dimensions, they respectively go like $y^{-2}$, $y^{-1}$ and $y^{0}$, up to logarithmic corrections. Further terms from the perturbative expansion of $F(\boldsymbol{\Delta})$ give contributions to $F(\boldsymbol{y})$ that vanish like $y$ or faster.

One may also derive the expansion (29) from the operator product expansion, without taking recourse to the transverse momentum representation (28). In the definition of collinear DPDs one has a product $\mathcal{O}_{2}(0,z_{2})\,\mathcal{O}_{1}(y,z_{1})$ of operators with $\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}$ but nonzero $\boldsymbol{y}$. This can be expanded around $\boldsymbol{y}=\boldsymbol{0}$ in terms of light-ray operators where all fields are at transverse position $\boldsymbol{0}$. These operators have twist 2, 3, 4 for the first, second and third term in (29), respectively.

The spitting contribution $F_{\text{spl,pt}}$ is given by graphs as in figures 4 and 5a and has already been discussed in the previous subsection. The term $F_{\text{tw3,pt}}$ originates from two types of graphs. The first type involves a single perturbative splitting and a quasipartonic collinear twist-three distribution as shown in figure 5b. The second type has two splittings as in figure 4 and a twist-three distribution with one “good” and one “bad” parton field. Given the helicity constraints discussed in section 2.1, collinear twist-three distributions in an unpolarised proton involve a quark and antiquark with opposite chirality (and possibly an extra gluon). As announced earlier, we discard twist-three terms in the following, since they are expected to become unimportant at small momentum fractions $x_{1},x_{2}$.

Finally, the term $F_{\text{int,pt}}$ contains contributions without any perturbative splitting; we hence refer to it as the “intrinsic” part of the DPD. It can be written as

[TABLE]

where $G$ is a quasipartonic collinear twist-four distribution and $C$ a perturbative splitting kernel, corresponding to graphs as in figure 6. The convolution $\otimes$ is in the longitudinal momentum fractions indicated by $\cdots$ (cf. figure 6a). The first term in (30) corresponds to the first term in (28) and is the only contribution that does not involve a hard splitting at all. The ellipsis denotes terms with non-quasipartonic twist-four distributions containing three or two parton fields, together with one or two parton splittings. While having the same power behaviour in $y$, one may expect that at small $x_{1},x_{2}$ these terms become less important than the terms with quasipartonic twist-four distributions, which should roughly grow as fast as the square of two parton densities.

We now take a closer look at the second term in (30). The kernel $C$ can be determined by computing both sides of (30) for a given graph. At order $\alpha_{s}$ only “non-diagonal” interactions, i.e. interactions connecting partons 1 and 2 as in figure 6a and b, contribute to $C$. The ladder graph in figure 6c is independent of $y$ and thus gives identical contributions to the matrix elements $F_{\text{int,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 (30). The l.h.s. evolves according to the homogeneous double DGLAP equation for DPDs, which describes “diagonal” interactions, either between the partons with final momentum fraction $x_{1}$ or between those with final momentum fraction $x_{2}$. By contrast, the evolution of $G(x_{1},x_{2},x_{2},x_{1};\mu)$ contains both diagonal and non-diagonal ladder interactions Bukhvostov:1985rn . The non-diagonal interactions in the evolution must thus be cancelled by the $\mu$ dependence of the term $C\otimes G$. At leading order in $\alpha_{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(\boldsymbol{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 (29) 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.

The following properties are in our opinion desirable for any theoretical setup describing double parton scattering.

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 $\boldsymbol{y}$ (or equivalently for all conjugate momenta $\boldsymbol{\Delta}$). 2. 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. 3.

For collinear factorisation, one would like to use as much as possible existing higher-order 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. 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 Diehl:2011yj ; Buffing:2016wip . 5. 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 Buffing:2016wip , adapting the well known procedure for single Drell-Yan production Collins:1984kg .

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 $\boldsymbol{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 Manohar:2012pe . However, DPDs and their products remain nonperturbative functions at large $\boldsymbol{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) for measured transverse momenta with $\Phi(y_{+}\nu)\,\Phi(y_{-}\nu)$ and the integrand in the collinear DPS formula (2) with $\Phi^{2}(y\nu)$. The function $\Phi(u)$ goes to $1$ for $u\gg 1$ and to [math] for $u\to 0$, and we can restrict ourselves to the case where $0\leq\Phi(u)\leq 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 Collins:2011zzd ). Consider a graph (or sum of graphs) $\Gamma$ that receives leading contributions from a set of loop momentum regions $R$. An approximation for $\Gamma$ is then given by

[TABLE]

In each term one integrates over all loop momenta. The operator $T_{R}$ applies approximations designed to work in momentum region $R$. Subtraction terms avoid double counting the contributions from smaller regions $R^{\prime}$ (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 Collins:2011zzd that $C_{R}\mskip 1.5mu\Gamma$ then provides a valid approximation in the region $R$ and in all smaller regions, and $\sum_{R}C_{R}\mskip 1.5mu\Gamma$ 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 $\Phi$ that depend on transverse distances $\boldsymbol{y}$ rather than transverse momenta. A collinear splitting region $R^{\prime}$ then corresponds to large $y$ and the corresponding hard region $R$ to small $y$, but we keep the ordering of regions $R^{\prime}<R$ from momentum space when implementing (31). 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}_{\pm}$ 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 $\Phi$ in the DPS cross section avoids an explicit modification of the definition of DPDs and thus respects criteria 3 and 4. In particular, the collinear DPDs $F(x_{i},\boldsymbol{y})$ in transverse position space follow the homogeneous DGLAP evolution equation (2). 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 $\nu$ to separate DPS from SPS and the twist-four contribution. In practical calculations one may take $\nu$ equal to the UV renormalisation scale $\mu$ in DPDs, but we find it useful to keep it separate in the general discussion. As a minimal requirement, $\nu$ must be of perturbative size, so that the double counting subtraction terms remove all contributions from nonperturbative regions in the hard kernels of the SPS and twist-four cross sections, making their calculation consistent. By construction, the $\nu$ dependence in the physical cross section cancels between DPS and the double counting subtraction terms, up to higher orders in $\alpha_{s}$ that are beyond the accuracy of the computation. We will come back to this issue in section 6.2.

Let us now take a closer look at the properties required for the function $\Phi$.

• •

The introduction of $\Phi$ in the DPS cross section must not spoil the correct description of the physics in the region where the DPD approximations work. Specifically, we require that the modifications introduced by $\Phi$ in that region should be power suppressed in the large scale. This requires that $\Phi(u)$ must approach $1$ sufficiently fast when $u\gg 1$.

Anticipating the momentum space analysis in section 4.3, we demand that for collinear factorisation the Fourier transformation of $1-\Phi^{2}$ exists in $2-2\epsilon$ dim (for positive and negative $\epsilon$). Therefore the integral \int du\,u^{1-2\epsilon}\bigl{[}1-\Phi^{2}(u)\bigr{]} must converge for $u\to\infty$. The corresponding criterion for TMD factorisation is obtained by replacing $\Phi^{2}\to\Phi$.

• •

$\Phi(u)$ at small $u$ must regulate the UV divergences of the naive DPD cross section in 2 dimensions. This means that $\int du\,u^{-1}\,\Phi(u)$ must be integrable at $u=0$ in TMD factorisation. For collinear factorisation we have the stronger requirement that $\int du\,u^{-3}\,\Phi^{2}(u)$ be integrable at $u=0$. These criteria are satisfied if for $u\to 0$

[TABLE]

with some $\delta>0$.

• •

To compute the double counting subtraction terms for SPS and the twist-four contribution, we must perform integrals over $\boldsymbol{y}$ in $2$ and $2-2\epsilon$ transverse dimensions, respectively, as we shall see in the next sections. We choose $\Phi$ such that the required integrals are known analytically. This is especially important for the twist-four contribution, where we must expand the result around $\epsilon=0$.

A suitable function for both collinear and TMD factorisation is the step function

[TABLE]

where $\gamma$ is Euler’s constant. This corresponds to a hard lower cutoff $y>b_{0}/\nu$ in the $y$ integration, where the constant $b_{0}\approx 1.12$ is taken to simplify certain analytical results. An alternative choice is

[TABLE]

Important integrals for these functions are given in table 1. Further possible choices for collinear factorisation are $\Phi^{2}(u)=1-\exp(-u^{p})$ or $\Phi^{2}(u)=u^{p}/(1+u^{p})$ with $p>2$; the corresponding integrals in the third and fourth row of table 1 can be performed and give Euler $\Gamma$ functions.

4.1 Leading order analysis: collinear factorisation

In this section we show how our formalism works in collinear factorisation, concentrating on the leading order in $\alpha_{s}$. Following the procedure of all-order factorisation proofs (see for instance Collins:2011zzd ), we rewrite individual Feynman graphs in a way consistent with the final factorisation formula. At the end of this procedure, the factors associated with long-distance physics in that formula, such as PDFs and DPDs, are expressed in terms of operator matrix elements and thus defined in a nonperturbative way. During our Feynman graph analysis we can separately consider splitting and intrinsic contributions to a DPD, $F_{\text{spl}}$ and $F_{\text{int}}$. Likewise, we can consider separate contributions $\sigma_{\text{1v1}}$, $\sigma_{\text{2v1}}$ and $\sigma_{\text{2v2}}$ to the DPS cross section, which respectively involve the combinations $F_{\text{spl}}\,F_{\text{spl}}$, $F_{\text{int}}\,F_{\text{spl}}+F_{\text{spl}}\,F_{\text{int}}$ and $F_{\text{int}}\,F_{\text{int}}$, corresponding to graphs a, b and c in figure 1. In the final factorisation formula, only the full DPDs $F$ and the full DPS cross section $\sigma_{\text{DPS}}$ will appear. For brevity, we write $\sigma$ instead of $d\sigma/(dx_{1}\mskip 1.5mudx_{2}\mskip 1.5mud\bar{x}_{1}\mskip 1.5mud\bar{x}_{2})$ here and in the following.

4.1.1 Squared box graph

Let us start with the squared box graph in figure 1a, integrated over all transverse momenta in the final state. As discussed in section 2, it receives its leading contribution from the region $q\sim Q$. The quark loops in the amplitude and its conjugate are then in the hard momentum region. This gives a part of the SPS cross section $\sigma_{\text{SPS}}$ and is computed as usual in terms of PDFs and the cross section for $gg$ annihilation into two gauge bosons. Other graphs that contribute to SPS, such as the one in figure 2a, have no overlap with DPS. In the present subsection we restrict $\sigma_{\text{SPS}}$ and other terms in the cross section to refer to the double box graph only.

The DPS region of the graph has near-collinear splittings in both the amplitude and its conjugate, which is only possible at $q\ll Q$. This gives the 1v1 contribution $\sigma_{\text{1v1}}$ to the DPS cross section, which is power suppressed compared with $\sigma_{\text{SPS}}$ because we integrate over the transverse boson momenta. Since the double box graph is integrated over all $q$ in $\sigma_{\text{SPS}}$, we subtract the contribution of the DPS region in order to prevent double counting. (At leading power accuracy this is not strictly necessary, but we will shortly see that it is useful). According to (31) this contribution is obtained by applying to the graph the approximations for DPS as well as those for SPS. This gives the DPS expression with the $g\to q\bar{q}$ splittings computed perturbatively, the incoming gluons treated as collinear on-shell partons, and quark masses in the loop set to zero. These approximations are just what goes into the perturbative splitting approximation $F_{\text{spl,pt}}$ of a DPD, which is given in (3.2). The subtraction term thus reads

[TABLE]

Note that in $\sigma_{\text{1v1,pt}}$ we use the short-distance approximation $F_{\text{spl,pt}}$ (even for large $\boldsymbol{y}$, where it is not valid), whereas in $\sigma_{\text{1v1}}$ we use the un-approximated splitting part $F_{\text{spl}}$ of each DPD. To ensure that $\sigma_{\text{1v1}}$ and $\sigma_{\text{1v1,pt}}$ receive contributions only from the DPS region, one should take $\nu\ll Q$.

To compute (35) we can use the known form (3.2) for the splitting DPDs, with the kernels collected in Diehl:2014vaa . With $\sigma_{\text{SPS}}$ being computed in fixed order perturbation theory, we make the same approximations in $\sigma_{\text{1v1,pt}}$, taking the renormalisation and factorisation scales in the splitting kernel and in the collinear PDF fixed. In $D=2$ dimensions the splitting DPDs then depend on $y$ like $1/y^{2}$. The $y$ integration in (35) is thus readily performed,

[TABLE]

This is proportional to $\nu^{2}$, which is expected since the unregulated integral has a quadratic divergence in $y$ and $\nu$ is the scale regulating this divergence. The relevant integral for different $\Phi(u)$ is obtained from the last row of table 1 by setting $\epsilon=0$.

We anticipate that the subtraction term at higher orders in $\alpha_{s}$ involves an integral as in (35), with splitting DPDs given by (3.2) and (27). Compared with (36) the integrand then has additional powers of $\log(y\mu)$. The integrals can be obtained from those in the last row of table 1 by Taylor expanding around $\epsilon=0$.

The complete contribution of the double box graph to the cross section is

[TABLE]

As is well known, an appropriate choice for the factorisation scale $\mu$ in collinear factorisation is $\mu\sim Q$. In the same spirit we take $\nu\sim Q$ rather than $\nu\ll Q$, although this extends the $y$ integrals in $\sigma_{\text{1v1}}$ and $\sigma_{\text{1v1,pt}}$ to values $y\sim 1/Q$ where the DPS approximation does not work. Both $\sigma_{\text{1v1}}$ and $\sigma_{\text{1v1,pt}}$ are then of the same order as $\sigma_{\text{SPS}}$ rather than power suppressed. Let us see that this still gives the correct approximation of the overall cross section.

Although $\sigma_{\text{SPS}}$ is naturally computed in momentum space, one can perform a Fourier transform w.r.t. the transverse momentum difference $\boldsymbol{\Delta}$ (defined as in figure 4 but for the graph of the cross section). Then $\sigma_{\text{SPS}}$ is a loop integral over the conjugate distance $\boldsymbol{y}$, just as the two other terms in (37).

We now show how the subtraction formalism works in this case. In the region of small $y\sim 1/Q$, one has $\sigma_{\text{1v1}}\approx\sigma_{\text{1v1,pt}}$, because the perturbative approximation of the DPD is designed to work at short distances. The dependence on the cutoff function $\Phi(y\nu)$ cancels between these two terms. One is therefore left with $\sigma_{\text{SPS}}$, which gives the appropriate description of the graph for $y\sim 1/Q$. We note that the DPDs in $\sigma_{\text{1v1,pt}}$ are computed at fixed order in $\alpha_{s}$, whereas in $\sigma_{\text{1v1}}$ they should be resummed using the double DGLAP equations. The cancellation between the two terms is therefore only up to higher orders in $\alpha_{s}$, which are beyond the accuracy of the result. We shall investigate this in detail in section 6.2.

Turning to the region of large $y\gg 1/Q$, one finds that $\sigma_{\text{SPS}}\approx\sigma_{\text{1v1,pt}}$, because precisely in that region the DPS approximation is designed to work for the graph. We are therefore left with the DPS term $\sigma_{\text{1v1}}$, which is exactly what we want to describe this region. The cutoff function $\Phi(y\nu)$ has no effect here, since we take $\nu\sim Q$ and require $\Phi(u)\approx 1$ for $u\gg 1$. The combination (37) thus gives the correct approximation for both small and large $y$.

4.1.2 2v1 graph

We now derive a description for the 2v1 graph in figure 1b. Integrated over transverse momenta in the final state, the graph receives leading contributions from small and large transverse momenta of the quarks on the splitting side (i.e. the top of the graph), corresponding to small or large transverse momenta of the produced gauge bosons. The region of small transverse momenta corresponds to DPS, whilst the contribution from large transverse momenta is correctly described by the twist-four mechanism introduced in section 2.

Whereas the SPS contribution in the previous subsection gives a finite result when computed as a perturbative two-loop graph, the twist-four contribution has a collinear divergence when integrated over the full region of phase space. To exhibit this divergence, we use dimensional regularisation and write

[TABLE]

As in the previous subsection, we restrict the meaning of $\sigma_{\text{tw4}}$ and other terms in the cross section to the contribution from a single graph. The labels $a_{1}$ and $a_{2}$ indicate parton species and polarisation, and $G$ is a collinear twist-four distribution. $H$ includes overall factors and the squared amplitude of the hard scattering process $q\bar{q}+g\to V_{1}V_{2}$. Note that in the hard scattering graph one has $\boldsymbol{q}_{1}=-\boldsymbol{q}_{2}=\boldsymbol{q}$. Momentum fractions are as shown in figure 7. The plus- and minus-momentum fractions of the produced bosons are respectively given by

[TABLE]

with $x_{i}$ and $\bar{x}_{i}$ defined in (4). For the momenta of internal lines in the hard scattering subprocess, minus and transverse components are fixed by the final state kinematics, so that there is only a loop integration over the plus-momentum fractions $v$ and $v^{\prime}$. Note that by virtue of (38) both $X_{i}$ and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{X}_{i}$ have an implicit dependence on $q^{2}$ at given $x_{i},\bar{x}_{i}$. An upper limit on $q^{2}$ follows from the requirements that $X_{1}+X_{2}\leq 1$ and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{X}_{1}+\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{X}_{2}\leq 1$.

The expression (4.1.2) includes a contribution from DPS region, characterised by $q_{i}\ll Q_{i}$ and $v,v^{\prime}\ll 1$. Let us approximate the integrand in that region. We first perform a Fierz transformation for the Dirac indices of the upper quark lines at the gauge boson vertices and retain only the leading Dirac structures. After this, the cross sections $\hat{\sigma}_{i}$ for $q\bar{q}$ annihilation appear in the expression. Furthermore we set $v=v^{\prime}=0$ in $G$ and extend the integrations over $v$ and $v^{\prime}$ to the full real axis (after which they are easily done using Cauchy’s theorem). Approximating $X_{i}\approx x_{i}$ and $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{X}_{i}\approx\bar{x}_{i}$ and collecting all factors, we obtain

[TABLE]

with the splitting kernel $T$ from (3.1). If we drop the restriction $q\ll Q_{i}$, then the expression in the second line of (4.1.2) becomes the perturbative splitting approximation $F_{\text{spl,pt}}(\bar{x}_{i},\boldsymbol{\Delta}=\boldsymbol{0})$ of the transverse-momentum integrated DPD. This DPD is evaluated at $\boldsymbol{\Delta}=\boldsymbol{0}$ in order to fulfil momentum conservation at the gauge boson vertices, because the transverse momenta of the right moving partons are set to zero in the overall hard scattering kernel $H$.

The $\boldsymbol{q}$ integral of (4.1.2) has an infrared divergence in $D=4$ dimensions, which must be cancelled by the subtraction term for the DPS region. That term is obtained by applying the DPS approximations to the graph, in addition to those for the twist-four region. This gives again the form (4.1.2), but with an unrestricted integration over $\boldsymbol{q}$, and a regulator of the corresponding ultraviolet divergence. To implement the regulator, we Fourier transform the collinear DPD, $F_{\text{spl,pt}}(\bar{x}_{i},\boldsymbol{\Delta}=\boldsymbol{0})=\int d^{D-2}\boldsymbol{y}\;F_{\text{spl,pt}}(\bar{x}_{i},\boldsymbol{y})$, and then multiply with $\Phi^{2}(y\nu)$ as prescribed by our formalism. This gives

[TABLE]

Alternatively, one can obtain the subtraction term by starting with the contribution $\sigma_{\text{2v1}}$ of graph 1b to the DPS cross section and applying the additional approximations adequate for the twist-four region. Let us translate the latter from transverse momentum (where the twist-four approximations are formulated) to transverse position (where our UV regulator for the DPS cross section is local). On the l.h.s. of the relation

[TABLE]

one neglects $\boldsymbol{\Delta}$ in the upper part of the graph, replacing $F_{\text{spl}}(\bar{x}_{i},\boldsymbol{\Delta})$ by its value at $\boldsymbol{\Delta}=\boldsymbol{0}$. The integral over $\boldsymbol{\Delta}$ of $F_{\text{int}}(x_{i},\boldsymbol{\Delta})$ gives the collinear twist-four distribution $G(x_{1},x_{2},x_{2},x_{1})$. In $\boldsymbol{y}$ space this corresponds to replacing $F_{\text{int}}(x_{i},\boldsymbol{y})$ by its value at $\boldsymbol{y}=\boldsymbol{0}$ while retaining the $\boldsymbol{y}$ dependence of $F_{\text{spl}}(\bar{x}_{i},\boldsymbol{y})$. Including the regulator function $\Phi^{2}(y\nu)$ under the $\boldsymbol{y}$ integral, one gets $F_{\text{int}}(x_{i},\boldsymbol{y}=\boldsymbol{0})\int d^{D-2}\boldsymbol{y}\;\Phi^{2}(y\nu)\,F_{\text{spl}}(\bar{x}_{i},\boldsymbol{y})$. The twist-four approximation implies computing the $g\to q\bar{q}$ splitting in perturbation theory, with massless quarks and the gluons taken collinear and on shell. This corresponds to replacing $F_{\text{spl}}$ with $F_{\text{spl,pt}}$. For the graph under discussion one has $F_{\text{int}}(x_{i},\boldsymbol{y}=\boldsymbol{0})=G(x_{1},x_{2},x_{2},x_{1})$ according to (30) and thus obtains (4.1.2). We see that, as in the case of the 1v1 graph, the subtraction term is obtained from the DPS cross section by replacing the DPDs with their small $y$ approximation at the appropriate order in $\alpha_{s}$. In section 5 we will find that this also holds at higher orders, where we have to take into account the second term in (30).

Let us finally perform the $\boldsymbol{y}$ integral in the subtraction term (4.1.2). Collecting all $\boldsymbol{y}$ and $\epsilon$ dependent factors in (3.2), we have

[TABLE]

where $\Omega_{2n}=2\pi^{n}/\Gamma(n)$ is the surface of a sphere in $2n$ dimensions. For brevity we do not display the momentum fraction in the splitting kernel $P$. The $\boldsymbol{y}$ integral is infrared divergent in $2$ dimensions and converges for $\epsilon<0$. Using the third row of table 1, we get in particular

[TABLE]

for $\Phi(u)=\Theta(u-b_{0})$, where $P^{\prime}=\partial P/\partial\epsilon$.

The contribution of the graph to the overall cross section is finally given by

[TABLE]

in analogy to (37). By the same argument as in the previous subsection, one finds that this combination reproduces the graph over the full range of $y$ values, with the second term cancelling against the third term for $y\sim 1/Q$ and against the first term for $y\gg 1/Q$.

As we have just seen, the product of $(\mu/\nu)^{2\epsilon}$ with the $1/\epsilon$ pole gives a logarithm $\log(\nu/\mu)$ in the subtraction term $\sigma_{\text{2v1,pt}}$. In turn, the $\boldsymbol{q}$ integral in the expression (4.1.2) of $\sigma_{\text{tw4}}$ gives $\log(Q/\mu)$, where $Q$ comes from the upper integration limit and the $\mu$ dependence from the regulated infrared divergence. In the combination $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ this gives $\log(Q/\nu)$. The $\nu$ dependence is cancelled by a $\log(\nu/\Lambda)$ from the $\boldsymbol{y}$ integral of the DPS cross section $\sigma_{\text{2v1}}$, where $\nu$ regulates the ultraviolet divergence. Combining all terms, one obtains $\log(Q/\Lambda)$ in (44). With the choice $\nu\sim Q$, the large part of the logarithm is contained only in the DPS term.

4.1.3 Combining contributions

We can now add all contributions to the physical cross section, including 1v1, 2v1 and 2v2 graphs in DPS. We include the same regulator $\Phi(y\nu)$ in all DPD contributions, so that the splitting and intrinsic contributions to the individual DPDs can be added up to distributions $F=F_{\text{spl}}+F_{\text{int}}$ that are defined by operator matrix elements as described in section 2.

2v2 graphs like the one in figure 1c are dominated by the DPS region, where transverse parton momenta are small, whilst regions with large transverse momenta are subleading. In impact parameter space this means that including the $\Phi$ regulator to suppress the region of small $y$ does not significantly change $\sigma_{\text{2v2}}$. We will see in section 4.3 that the effect of $\Phi$ in the 2v2 term is at the level of power corrections.

The master formula for the cross section with integrated transverse momenta is then

[TABLE]

where $\sigma_{\text{DPS}}=\sigma_{\text{1v1}}+\sigma_{\text{2v1}}+\sigma_{\text{2v2}}$ includes all contributions to the DPS cross section and involves the full DPDs $F$. From the discussion in the previous two subsections it follows that the $\nu$ dependence of $\sigma_{\text{DPS}}$ cancels against the one of $\sigma_{\text{1v1,pt}}+\sigma_{\text{2v1,pt}}$, up to terms of higher order in $\alpha_{s}$ that are beyond the accuracy of the calculation. As discussed in section 2.1 we neglect contributions involving twist-three functions in (45).

The terms in (45) are now meant to include all contributing graphs, including graphs without a double counting issue. The hard scattering cross sections needed for the DPS term are available for many final states, and $\sigma_{\text{1v1,pt}}$ can be readily obtained from existing results at leading order in $\alpha_{s}$. For double gauge boson production, the SPS cross section has been computed at high perturbative order, including the double box graphs as well as the first radiative corrections to them Caola:2016trd . On the other hand, to the best of our knowledge, no calculations of $\sigma_{\text{tw4}}$ have been performed so far. To leading logarithmic accuracy, one can however omit the terms $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ in (45), because with the choice $\nu\sim Q$ the large DPS logarithm $\log(Q/\Lambda)$ generated by 2v1 graphs is entirely contained in $\sigma_{\text{DPS}}$, as discussed in the previous subsection.

4.2 Leading order analysis: TMD factorisation

We now show how our formalism works for TMD factorisation. As explained in section 2.1, the cross section for $q_{1},q_{2}\ll Q$ receives leading-power contributions from DPS, SPS and from their interference. The 2v1 contribution is power suppressed in this case and will not be discussed further in this subsection. The discussion of leading momentum regions is a bit more involved now, because the two quark loops in the double box graph can be collinear or hard independently of each other.

We start with the SPS/DPS interference, for which a graph is shown in figure 2b. Its contribution to the $V_{1}V_{2}$ production cross section has the form

[TABLE]

for DPS in the amplitude or in its complex conjugate. Here $D$ are the twist-three TMDs introduced in section 3.1 and the proportionality is up to kinematic and numerical factors. We have different hard scattering amplitudes $H$, with $\alpha_{0}$ and $\beta_{0}$ being gluons and the remaining partons quarks or antiquarks . The regions of small $y_{+}$ or $y_{-}$ are regulated by $\Phi(y_{+}\nu)$ and $\Phi(y_{-}\nu)$.

Notice that the graph in figure 2b has a leading contribution when the $g\to q\bar{q}$ splittings in the conjugate amplitude become collinear, in which case we are back in the DPS region. From the first expression in (4.2) one should therefore subtract

[TABLE]

with $F_{y_{-}\to 0}$ given in (19) and represented in figure 5b. From the second expression in (4.2) one should subtract an analogous term with $y_{+}\to 0$.

According to (18), the computation of the subtraction terms involves the integrals

[TABLE]

which are given in the appendix. We see in table 1 that for $n=0$ both our choices for $\Phi$ give a logarithmic behaviour $\log(\nu/q)$ when $\nu\gg q$, as expected for a regulated logarithmically divergent integral.

Let us now turn to the double box graph. It has leading contributions when the loops are both hard (SPS region) or both collinear (DPS region), or when one loop is hard and the other one collinear (region of SPS/DPS interference). The terms to be added to the SPS cross section for removing the contribution from the SPS/DPS interference regions are

[TABLE]

The first terms in the square brackets are obtained from (4.2) by taking the perturbative splitting form for all distributions, the one for $D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})$ being given in (21). The second terms are obtained from (4.2) and its analogue for $y_{+}\to 0$ by taking both splitting vertices perturbative, using the DPD given in (3.1). Notice that the subtraction terms (4.2) contain subtractions themselves, showing the recursive nature of the construction in (31). Finally, one should remove the contribution from the DPS region by adding

[TABLE]

where again both splitting vertices in each DPD are taken as perturbative. Note that in both (4.2) and (50) all proton matrix elements are expressed in terms of ordinary twist-two TMDs, as is the case for the SPS cross section.

Adding up contributions and double counting subtractions, we finally obtain

[TABLE]

where we have again omitted differentials $dx_{1}\,dx_{2}\,\ldots$ for brevity. The DPS cross section $\sigma_{\text{DPS}}$ contains contributions with and without splitting in each of the DPDs, which are defined by operator matrix elements as in the case of collinear factorisation.

Let us now discuss the large logarithms appearing in the different terms. The hard scattering kernel $H_{\alpha_{0}\beta_{0}}$ in (4.2) has a DPS logarithm $\log(Q/q)$, and correspondingly there is a $\log^{2}(Q/q)$ in the unsubtracted SPS cross section. We take $\nu\sim Q$, so that these logarithms are fully removed by the double counting subtractions. They reappear in the logarithmic integrals over $\boldsymbol{y}_{+}$ and $\boldsymbol{y}_{-}$ in the DPS cross section and in the SPS/DPS interference terms (4.2). In the full cross section (4.2) one then has a squared logarithm $\log^{2}(Q/q)$ in the DPS term, and a single $\log(Q/q)$ in the subtracted SPS/DPS interference in the first line, whereas the subtracted SPS contribution in the last line has no leading logarithm. If $q$ is in the nonperturbative region, one has $\log(Q/\Lambda)$ instead of $\log(Q/q)$.

At a practical level, the computation of the SPS/DPS interference requires twist-three TMDs as input, and the DPS term requires transverse-momentum dependent DPDs. In general, the modelling of these functions is very difficult and largely unconstrained. The situation is much better when $q_{1},q_{2}\gg\Lambda$: in this case one can compute these functions in terms of collinear ones Buffing:2016wip .

On the perturbative side, the calculation of the hard scattering cross sections in the DPS term is easiest, while the SPS and interference terms require the full computation of box graphs. If one is satisfied with keeping only the leading double logarithm in $Q/q$, then it is sufficient to retain only the DPS contribution with $\nu\sim Q$.

4.3 Subtraction formalism in momentum space

In Collins:2011zzd it was proven that the subtraction formalism leading to (31) correctly approximates a graph up to power suppressed terms, provided that each approximant $T_{R}$ correctly reproduces the graph in its design region $R$, up to power suppressed terms. This proof works with regions in momentum space. However, our UV regulator $\Phi$ for DPS singularities is multiplicative in $\boldsymbol{y}$ space, which gives convolutions in transverse momenta that do not appear in the original Feynman graphs and hence are not part of the proof in Collins:2011zzd . We now show that this does not cause problems, so that our setup gives adequate approximations when transformed to momentum space, with corrections that are subleading in powers of $\Lambda/\nu$. With $\nu\sim Q$ these corrections are of the same order as other approximations in the factorisation formula.

We begin with the case of collinear factorisation, using dimensional regularisation for UV and IR divergences as usual. As mentioned earlier, we choose $\Phi$ such that the Fourier transform

[TABLE]

is finite for all values of $\Delta/\nu$. In the following we show that terms that depend on $\Psi$ in the cross section are power suppressed and can hence be discarded when establishing the validity of the subtraction formalism.

Let us first discuss $\sigma_{\text{2v2}}$, which is proportional to

[TABLE]

The product $\hat{\sigma}_{1}\mskip 1.5mu\hat{\sigma}_{2}$ of hard scattering cross sections does not affect our argument and is hence omitted. To keep the notation simple, we only display transverse-momentum arguments in this subsection. For simplicity we perform the power counting for $D=4$ dimensions; changes in $D=4-2\epsilon$ are by fractional powers and do not alter our conclusions. We then have

[TABLE]

for both $\Delta\sim\Lambda$ and $\Delta\sim\nu$, which are the two scales present in (4.3).555In the ultraviolet region $\Delta\gg\nu$ there is a further suppression of $\Psi(\boldsymbol{\Delta})$ compared with $1/\nu^{2}$, the exact form of which depends on the function $\Phi$. The large $\Delta$ behaviour of $F_{\text{int}}$ follows from dimensional analysis ($\Lambda^{2}$ in the numerator comes from a collinear twist-four matrix element and $1/\Delta^{2}$ from the hard splitting kernel in the Fourier transformed version of (30)). The $\Psi$ independent term in (4.3) receives its leading contribution of order $\Lambda^{2}$ from the region $\Delta\sim\Lambda$, which is the design region of the DPS approximation. By contrast, the $\Psi$ dependent term behaves like $\Lambda^{4}/\nu^{2}$, with contributions from both $\Delta\sim\Lambda$ and $\Delta\sim\nu$. We thus see that our UV regulator (which is not necessary in $\sigma_{\text{2v2}}$ as noted earlier) only gives changes suppressed by $\Lambda^{2}/\nu^{2}$. For $\nu\sim Q$ this does not degrade the overall accuracy of the calculation.

Next we consider the contribution $\sigma_{\text{2v1}}-\sigma_{\text{2v1,pt}}+\sigma_{\text{tw4}}$ from 2v1 graphs. Here we have a subtraction term, which is obtained by replacing the DPDs in $\boldsymbol{y}$ space with their perturbative approximations. In momentum space this reads

[TABLE]

The second line is obtained by replacing $F$ with $F_{\text{int}}$ in the first term and with $F_{\text{int,pt}}$ in the second term on the r.h.s. of (28) and then taking the Fourier transform with respect to $\boldsymbol{y}$. We thus find that $\sigma_{\text{2v1}}-\sigma_{\text{2v1,pt}}$ is proportional to

[TABLE]

For power counting we use (54) and additionally

[TABLE]

which follows from dimensional analysis of the splitting graph. The perturbative approximations $F_{\text{int,pt}}$ and $F_{\text{spl,pt}}$ scale like their un-approximated counterparts. The $\Psi$ independent terms in (4.3) are then found to be of order $\Lambda^{2}$, which after multiplication with $\hat{\sigma}_{1}\mskip 1.5mu\hat{\sigma}_{2}$ gives the same power behaviour in the cross section as $\sigma_{\text{tw4}}$. To analyse the $\Psi$ dependent terms, we rearrange them as follows:

[TABLE]

The first line is suppressed by the integration phase space for $\Delta\sim\Lambda$, by $\Psi(\boldsymbol{\Delta}-\boldsymbol{\Delta}^{\prime})-\Psi(\boldsymbol{\Delta})$ for $\Delta\sim\nu$ and $\Delta^{\prime}\sim\Lambda$, and by $F_{\text{int}}(-\boldsymbol{\Delta}^{\prime})-F_{\text{int,pt}}(-\boldsymbol{\Delta}^{\prime})$ for $\Delta,\Delta^{\prime}\sim\nu$. The second line is again suppressed by the integration phase space for $\Delta\sim\Lambda$, and by $F_{\text{spl}}(\boldsymbol{\Delta})-F_{\text{spl,pt}}(\boldsymbol{\Delta})$ for $\Delta\sim\nu$. The sum of all $\Psi$ dependent terms is hence power suppressed compared with the $\Psi$ independent terms. One readily finds that the relative suppression is by $\Lambda^{2}/\nu^{2}$, as it was in the 2v2 case. Notice that this suppression is obtained for the combination $\sigma_{\text{2v1}}-\sigma_{\text{2v1,pt}}$ but not for the two terms separately.

Finally, we discuss the contribution $\sigma_{\text{1v1}}-\sigma_{\text{1v1,pt}}+\sigma_{\text{SPS}}$ from 1v1 graphs. The combination $\sigma_{\text{1v1}}-\sigma_{\text{1v1,pt}}$ is proportional to

[TABLE]

One finds that both the $\Psi$ dependent and the $\Psi$ independent terms on the r.h.s. scale like $\Lambda^{2}$. Multiplied with $\hat{\sigma}_{1}\mskip 1.5mu\hat{\sigma}_{2}$ this is power suppressed by $\Lambda^{2}/Q^{2}$ compared with $\sigma_{\text{SPS}}$. Notice that the SPS approximation itself has relative corrections of order $\Lambda^{2}/Q^{2}$. This means that the contribution in (4.3) is beyond the accuracy of the calculation of 1v1 graphs. The fact that the $\Psi$ dependent part is of similar size as the $\Psi$ independent one is hence of no concern. We recall from section 2.1 that the rationale to evaluate the 2v2 and 2v1 terms, which are also suppressed by $\Lambda^{2}/Q^{2}$, is that they may be enhanced by other factors.

The case of TMD factorisation can be discussed along the same lines. In analogy to (52) we define $\Psi$ as the Fourier transform of $[\mskip 1.5mu1-\Phi(y\nu)]$ instead of $[\mskip 1.5mu1-\Phi^{2}(y\nu)]$. Let us discuss the SPS/DPS interference graph in figure 2b. Its contribution to $\sigma_{\text{DPS/SPS}}$ in (4.2) is proportional to

[TABLE]

where on the r.h.s. we have twist-three TMDs in momentum representation. An analogous relation holds for the contribution to $\sigma_{\text{DPS/SPS},\,y_{+}\to 0}$, which according to (4.2) is to be subtracted from $\sigma_{\text{DPS/SPS}}$ in the overall cross section. In this case one should take the perturbative splitting approximation $D_{\text{spl,pt}}$ of the distributions, which in $\boldsymbol{y}_{+}$ representation is given by (21). For power counting we take $q_{1},q_{2}\sim\Lambda$ and use

[TABLE]

where $k_{-}$ and $K$ can be of order $\Lambda$ or $\nu$. This can be derived by a general analysis of graphs as in section 5.2 of Diehl:2011yj ; the $1/k_{-}$ dependence is also directly obtained by Fourier transforming (21). $D_{\text{spl,pt}}$ behaves like $D$, and of course we have $\Psi\sim 1/\nu^{2}$ as before.

We see that for $k_{-}\sim K\sim\Lambda$ the $\Psi$ independent term in (4.3) is of order $1/\Lambda^{2}$, whilst the $\Psi$ dependent one is of order $1/\nu^{2}$. The only region in which the $\Psi$ dependent term gives a leading contribution is when $k_{-}\sim k^{\prime}_{-}\sim\nu$ and $K\sim\Lambda$. However, in this region one has $D\approx D_{\text{spl,pt}}$, so that its contribution is suppressed in the combination $\sigma_{\smash{\text{DPS/SPS}}}-\sigma_{\smash{\text{DPS/SPS},\,y_{+}\to 0}}$. We thus find that in the overall cross section, $\Psi$ dependent terms are power suppressed, as in the case of collinear factorisation. The same holds of course for $\sigma_{\smash{\text{SPS/DPS}}}-\sigma_{\smash{\text{SPS/DPS},\,y_{-}\to 0}}$. Finally, a corresponding analysis (involving the momenta $\boldsymbol{k}_{+}$, $\boldsymbol{k}_{-}$ and $\boldsymbol{K}$) can be given for the combination $\sigma_{\text{DPS}}-\sigma_{\text{DPS},\,y_{-}\to 0}-\sigma_{\text{DPS},\,y_{+}\to 0}+\sigma_{\text{DPS},\,y_{\pm}\to 0}$, which covers all terms in (4.2) that depend on our UV regulator.

5 Subtraction terms at higher orders

Our construction of the double counting subtractions for DPS regions is not limited to the leading-order graphs discussed so far. In the present section we illustrate how the formalism works at highers in $\alpha_{s}$, taking as an example the 2v1 mechanism in collinear factorisation. In the master formula (31) we encounter nested subtractions in this case. We show for selected graphs how the formalism gives the different contributions to the cross section, namely $\sigma_{\text{tw4}}$, $\sigma_{\text{2v1}}$ and the subtraction term $\sigma_{\text{2v1,pt}}$. In particular we find that the latter can be obtained from $\sigma_{\text{2v1}}$ by replacing the DPDs with their short-distance expansions $F_{\text{int,pt}}$ and $F_{\text{spl,pt}}$, as found at leading order in section 4.1.2.

Figure 8: LO and NLO graphs for the 2v1 mechanism, with the blobs for the hadronic matrix elements omitted as in figure 7.

We focus here on the general structure and use a schematic notation, in particular suppressing indices and arguments that are not essential for the discussion.

Graph 5a.

To introduce this notation, let us briefly review the leading-order graph in figure 5a, which was discussed in detail in section 4.1.2. It has two leading momentum regions, which we denote by $C$ and $H$. In region $C$ the gauge bosons have small transverse momenta and the $g\to q\bar{q}$ splitting is near collinear. This is the region of DPS, and its approximated contribution to the cross section reads

[TABLE]

Here $\hat{\sigma}^{(0)}_{i}$ is the tree-level cross section for $q\bar{q}\to V_{i}$, and we write $\int_{y}$ as a shorthand for $\int d^{D-2}\boldsymbol{y}$. The contribution of this graph to the splitting DPD at the top is denoted by $F^{(0)}_{\text{spl}}$, and its contribution to the intrinsic DPD at the bottom by $F^{(0)}_{\text{int}}$. Here and in the rest of the section, superscripts $(0),(1),(2)$ indicate the power of $\alpha_{s}$ in the considered quantity. For simplicity we suppress longitudinal momentum arguments throughout.

In region $H$, the transverse momenta of the gauge bosons are large and thus the $q\bar{q}$ pair into which the gluon splits is far off shell on both sides of the final state cut. This is described by the twist-four contribution and requires a subtraction for the smaller region $C$. The subtracted contribution is

[TABLE]

Here $\int_{q}$ is shorthand for $\int d^{D-2}\boldsymbol{q}$, the tree-level hard scattering cross section for $q\bar{q}+g\to V_{1}V_{2}$ is denoted by $H^{(1)}$, whilst $f^{(0)}$ and $G^{(0)}$ respectively are the zeroth order contributions to the gluon distribution at the top and the twist-four $q\bar{q}$ distribution at the bottom of the graph. As seen in (4.1.2), there is a convolution integral over longitudinal momentum fractions in $G$, denoted here by $\otimes$. In the subtraction term, $V^{(1)}$ represents the lowest-order perturbative $g\to q\bar{q}$ splitting kernel, including all terms in (3.2) except for the single-parton density $f_{a_{0}}$. In this term, the momentum fractions of $G^{(0)}$ are fixed by external kinematics, so that the convolution $\otimes$ is absent.

The perturbative expressions of the DPDs for small $y$ are

[TABLE]

at leading order, so that the combined contribution of the $C$ and $H$ regions can be written as

[TABLE]

As already mentioned, the double counting subtraction in the second term is simply obtained from the DPS term by replacing each DPD with its small-$y$ expansion at the appropriate order in $\alpha_{s}$.

Graph 5b.

We now analyse the NLO graph in figure 5b. The $q\bar{q}$ splitting can be hard or collinear, and likewise the additional gluon with momentum $\ell$. The leading momentum regions thus are $CC$, $HC$ and $HH$, where the first letter refers to the $q\bar{q}$ splitting and the second to the additional gluon. There is also a region where the gluon is hard whereas the $q\bar{q}$ splitting to the right of the cut is collinear. The splitting at the left of the cut is then hard by momentum conservation. One finds that this region is not leading: the loss of phase space for having the gauge bosons at low transverse momenta is stronger than the gain from having one but not both splitting vertices collinear. Let us discuss the leading regions in turn.

In the $CC$ region, the extra gluon is part of the DPD for the right moving particles, and one simply has

[TABLE]

where the superscript $(1)$ indicates a contribution to $F_{\text{int}}$ at order $\alpha_{s}$. Likewise, in the $HC$ region, the extra gluon becomes part of the twist-four distribution $G$, so that one has

[TABLE]

Note that the integral over the gluon momentum in the graph giving $G^{(1)}_{\text{ren}\rule{0.0pt}{3.31529pt}}$ contains an ultraviolet divergence. It is understood that $G^{(1)}_{\text{ren}\rule{0.0pt}{3.31529pt}}$ is renormalised, i.e. includes the appropriate UV counterterm, which is indicated by the subscript “ren”.

For the $HH$ region, the recursive character of the subtraction formalism becomes evident. We have a subtraction for the next smallest region $HC$, which itself contains a subtraction for the region $CC$. In addition, a $CC$ subtraction has to be made for the overall graph. This gives

[TABLE]

The term in the first line is for the hard region, with $H^{(2)}$ being a one-loop contribution to the cross section for $q\bar{q}+g\to V_{1}V_{2}$. In the second line, the term for the $HC$ region (including its subtraction for the $CC$ region) involves a non-diagonal splitting kernel $\widetilde{K}^{(1)}_{\text{ren}}(0)$ for the exchange of a gluon between the two parton lines with momentum fractions $x_{1}$ and $x_{2}$, which are taken at relative transverse position $y=0$ because this is part of the approximation for the $HC$ region. In the last line, which is the overall subtraction for the $CC$ region, the same kernel is to be taken at relative distance $y$, because in this case the $q\bar{q}$ splitting is not assumed to be harder than the extra gluon emission. One therefore has to take a convolution in transverse momenta, which turns into a product of two $y$ dependent factors $V^{(1)}(y)$ and $\widetilde{K}^{(1)}(y)$ after Fourier transformation.

Notice that in dimensional regularisation the limit $y\to 0$ in $\widetilde{K}^{(1)}(y)$ is not smooth. The graph for $\widetilde{K}^{(1)}_{\text{ren}}(0)$ contains a scaleless momentum integral, which is zero. After renormalisation, i.e. subtraction of the ultraviolet counterterm, one is left with an infrared divergence.666The UV divergent part of $\widetilde{K}^{(1)}_{\text{ren}}(0)$ gives DGLAP splitting kernels generalised to non-diagonal momentum fractions, known from the evolution of twist-four distributions Bukhvostov:1985rn and of GPDs Belitsky:2005qn . By contrast, $\widetilde{K}^{(1)}(y)$ at finite $y$ does not have an ultraviolet divergence (and hence no counterterm), but it does have the same infrared divergence as $\widetilde{K}^{(1)}_{\text{ren}}(0)$.

Adding all contributions and using (64) we obtain

[TABLE]

with

[TABLE]

In the first line of (5) we have an NLO contribution to the twist-four term, with the NLO kernel $H^{(2)}$ having a subtraction for the region where the gluon becomes collinear, as is standard in collinear factorisation. In the DPS subtraction term we have a contribution to the small-$y$ expansion of $F_{\text{int}}$ at NLO, with the three terms in (71) corresponding respectively to the three terms in the general expression (28). The infrared divergence of the kernel $K$ drops out here, as it must. We see that the DPS subtraction term in (5) has the same form as in (5), with $F_{\text{int,pt}}^{(1)}$ instead of $F_{\text{int,pt}}^{(0)}$.

In (5) all ultraviolet and infrared subtractions have been carried out, so that one can set $D=4$ in dimensional regularisation. The quantities with subscript “ren” then depend on the associated scale $\mu$. In the first line, this dependence cancels between $G^{(1)}_{\text{ren}\rule{0.0pt}{3.31529pt}}$ (the graph under discussion contributes to the scale evolution of $G$) and the subtracted hard scattering kernel in the second term. In (71), the $\mu$ dependence cancels between $G^{(1)}_{\text{ren}\rule{0.0pt}{3.31529pt}}$ and the subtracted kernel of the small-$y$ expansion, i.e. graph 5b does not contribute to the evolution of $F_{\text{int}}$.

Graph 5c.

We next consider the graph in figure 5c. The contribution from the region $CC$ is

[TABLE]

where in contrast to graph 5b the one-loop contribution to $F_{\text{int}}$ now has an ultraviolet divergence (since the two partons connected by the gluon are at the same transverse position) and needs renormalisation. The contribution from region $HC$ is analogous to (68). Unlike for graph 5b, a leading contribution is now also obtained for the momentum region $CH$, where the $g\to q\bar{q}$ splittings are collinear, whilst the gluon carries a large transverse momentum (balanced by the vector boson momentum $\boldsymbol{q}_{2}$). Including the subtraction for the region $CC$, we have

[TABLE]

where $\sigma^{(1)}_{2}$ is the cross section for $q\bar{q}\to V_{2}+g$ and $K^{(1)}_{\text{ren}}$ is a diagonal splitting kernel including the necessary UV counterterm. Its divergent part involves the contribution of this graph to the usual DGLAP splitting kernel $P_{qq}$. Finally, the region $HH$ gives

[TABLE]

where the last three lines are the subtraction terms for the regions $HC$, $CH$, $CC$, respectively. Adding all contributions, we obtain

[TABLE]

where in the second line we used (64), (65) and

[TABLE]

The latter relation holds because the exchanged gluon in graph 5c gives exactly the same contribution to $F_{\text{int}}(y)$ and to $G$ and hence cancels between the second and third term in the general small-$y$ expansion (28) of $F_{\text{int}}(y)$. In each line of (5) we have hard scattering cross sections subtracted for the collinear region, as is characteristic for collinear factorisation. The $\mu$ dependence cancels between the subtracted cross sections and the renormalised distributions.

Graph 5d.

We now discuss figure 5d, where we have an additional gluon in the upper part of the graph. Its leading momentum regions are $CC$, $HC$ and $HH$. As in the case of graph 5b, there is no leading region $CH$. For the $CC$ region we simply have

[TABLE]

In the region $HC$, the extra gluon in the graph is collinear and couples to an off-shell quark, since the $g\to q\bar{q}$ splitting is hard. Among the approximations to be made in this case is the Grammer-Yennie approximation (see e.g. Diehl:2011yj ; Diehl:2015bca ). After summing over an appropriate set of graphs, one can apply a Ward identity, after which the gluon couples to an eikonal line, as shown in figure 9. We thus can write

[TABLE]

where $\Gamma^{e+}$ denotes the sum of graph e and of two graphs analogous to graphs d and e, with the gluon coupling to the rightmost quark line. $f^{(1)}_{\text{ren}}$ is an NLO contribution to the gluon density, involving a gluon coupling to an eikonal line. The integral over the gluon transverse momentum diverges in the ultraviolet and hence requires renormalisation.

For the $HH$ region we have

[TABLE]

In the subtraction term for the region $HC$, we have an NLO contribution $K^{(1)}_{\text{ren}}$ to the $gg$ splitting kernel with a gluon coupling to an eikonal line, in analogy to $f^{(1)}_{\text{ren}}$ above. Again, this requires summation over several graphs as specified on the l.h.s. In the last line of (5) we have the subtraction term for the region $CC$, in which the gluon and the splitting are regarded as collinear, so that no eikonal approximation is made and one obtains an NLO contribution $V^{(2)}$ to the $g\to q\bar{q}$ splitting kernel. Adding all contributions, we obtain

[TABLE]

where

[TABLE]

is a contribution to the perturbative splitting DPD at NLO. The scale dependence of $f^{(1)}_{\text{ren}}$ and $K^{(1)}_{\text{ren}}$ cancels in (81), corresponding to the fact that graph 5d and the other ones required for applying the Ward identity do not contribute to the scale evolution of the DPD.

Other NLO graphs.

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

6 Double DGLAP evolution

Collinear DPDs $F(x_{i},\boldsymbol{y};\mu_{i})$ evolve with $\mu_{i}$ according to the double DGLAP equations. In the present section, we discuss several aspects of this scale dependence.

6.1 Resummation of large logarithms

As we have seen in section 4.1.2, the 2v1 graph in figure 5a gives rise to a large logarithm $\log(Q/\Lambda)$ when the cross section is integrated over the transverse boson momenta. This logarithm originates from the $g\to q\bar{q}$ splitting vertices in the kinematic region of double parton scattering and is correctly reproduced in the DPS cross section, provided that we make a suitable choice of scales.

Radiative corrections to figure 5a give rise to additional large logarithms of DGLAP type. One can expect that these logarithms are correctly resummed to all orders by the DGLAP evolution of the DPDs, provided again that one makes suitable scale choices. In this section we show for explicit examples that this is indeed the case. We will allow for different hard scales $Q_{1},Q_{2}$ and set the scales as

[TABLE]

where we define

[TABLE]

Of course, our conclusions regarding large logarithms will not change if (82) is modified by factors of order unity. For definiteness, we will always take a rigid cutoff $\Phi(u)=\Theta(u-b_{0})$ in $y\nu$. For later use we define the scale

[TABLE]

DGLAP logarithms are generated by loop integrals in regions of strongly ordered transverse momenta. We will compare the logarithms obtained directly from transverse momentum integrals with the ones obtained from evolved DPDs in $y$ space.

To set the scene, we start with the lowest-order 2v1 graph in figure 5a. In the DPS region, i.e. when the scale of the $g\to q\bar{q}$ splittings is much smaller than $Q_{1}$ and $Q_{2}$, the behaviour of a splitting vertex is easily recovered from the computation of the splitting DPD. As can be seen from (3.1), a vertex with quark momentum $k$ gives a factor $\boldsymbol{k}^{j}/\boldsymbol{k}^{2}$, where $j$ is a transverse index. The transverse momentum of the splitting vertices must flow through the gauge boson lines. Thus, the two splitting vertices in figure 5a give a combined $1/q_{1}^{2}$ behaviour for $q_{1}\ll Q_{\text{min}}$, as is evident in (4.1.2) after angular integration. The integrated cross section is then proportional to

[TABLE]

Since the $q_{1}$ integral is logarithmic, we must include an infrared cutoff $\Lambda$, which limits the integration to the perturbative region where the preceding analysis is valid. We take here the usual procedure for extracting the leading behaviour of a logarithmic integral $\int_{a}^{b}dk\;\Gamma(k)$, approximating the integrand $\Gamma(k)$ for $a\ll k\ll b$ while integrating over the full range $a\leq k\leq b$.

Let us now see how the logarithm (85) arises in the DPS cross section formula. According to the perturbative expressions of the DPDs, $F_{\text{spl,pt}}\sim 1/y^{2}$ and $F_{\text{int,pt}}$ is $y$ independent at leading order in $\alpha_{s}$, so that we have

[TABLE]

With our choice for the infrared cutoff on $y$, the expressions (85) and (86) agree including their overall normalisation.

Now consider the NLO graph in figure 10a. We route the gluon momentum $\ell$ back through the boson line with momentum $q_{2}$, so that $\ell$ flows through exactly two quark propagators. In the strongly ordered region

[TABLE]

the graph has the same DPS logarithm as the leading-order result, whilst developing an additional logarithm from the integral $\int d\ell^{2}/\ell^{2}$, as is familiar from similar gluon emission graphs in hard scattering processes. In the case $Q_{2}<Q_{1}$ this gives

[TABLE]

whereas for $Q_{1}<Q_{2}$ we have

[TABLE]

In both cases, we have an additional logarithm $\log(Q_{\text{min}}^{2}/\Lambda^{2})$ compared with the lowest-order result. If $Q_{1}\ll Q_{2}$, then the logarithm of the ratio $Q_{2}/Q_{1}$ is also large.

Let us compare these results with the 2v1 part of the DPS cross section formula when we take into account evolution of the splitting DPD. Focusing on the dependence on $y$ and $\mu$, we can write

[TABLE]

at leading order, where the choice of scale $\mu_{y}$ ensures that there are no large logarithms of $y\mu$ in the higher-order corrections. This should be evolved to $F_{\text{spl,pt}}(x_{i},\boldsymbol{y};\mu_{1},\mu_{2})$ in the cross section. In (90) there is a logarithmic dependence on $y$ from the running of $\alpha_{s}$ and from the scale dependence of the gluon density. Graph 10a does not contribute to either of these, so that the corresponding dependence should not be taken into account here. Rather, the DGLAP logarithm of graph 10a corresponds to the evolution of (90) from $\mu_{y}$ to $\mu_{2}$, which gives a factor proportional to $\alpha_{s}\log(\mu_{2}^{2}/\mu_{y}^{2})$. The relevant $y$ integral in the DPD cross section thus is

[TABLE]

where we have used the scale choice in (82). This agrees with both (88) and (89).

Note the ease with which (91) is obtained in $y$ space, compared with the detailed considerations leading to (88) and (89). We see that with evolved DPDs in $\sigma_{\text{DPS}}$ one includes higher order logarithms while using tree-level partonic cross sections. Treating the 2v1 graphs fully within the collinear twist-four formalism would require the explicit inclusion of the higher-order graphs.

Whereas in the previous example, the gluon is emitted close to the splitting vertices, the opposite is true for the graphs in figure 10b and c. For figure 10b the region generating two large logarithms is

[TABLE]

which gives

[TABLE]

To see how these logarithms arise from the evolution of $F_{\text{int,pt}}$, we use the short-distance expansion

[TABLE]

at leading order, which we evaluate at scale $\mu_{y}$ to avoid large logarithms of $y\mu$ in the higher-order corrections, as we did for $F_{\text{spl,pt}}$ in (90). Graph 10b contributes to the scale evolution of $G$, where at leading order it gives $\alpha_{s}\log(\mu_{y}^{2}/\Lambda^{2})$ when evolution is started from the generic soft scale $\Lambda^{2}$. The graph also contributes to the evolution of $F_{\text{int,pt}}$ from $\mu_{y}$ to $\mu_{2}$, where it gives $\alpha_{s}\log(\mu_{2}^{2}/\mu_{y}^{2})$. The two effects add up to a factor $\alpha_{s}\log(\mu_{2}/\Lambda^{2})$. In the DPS cross section we thus have

[TABLE]

in agreement with (93).

The graph in figure 10c also has a region where DGLAP logarithms are generated, as was already discussed in Gaunt:2012dd . We route the loop momentum $\ell$ through two quark lines at one of the two $g\to q\bar{q}$ splitting vertices. In the strongly ordered region

[TABLE]

the two lower fermion lines to which the gluon couples have virtualities of order $\ell^{2}$, whereas both splitting vertices have virtualities of order $q_{1}^{2}\approx q_{2}^{2}$. This gives 777Note that in equation (3.13) of Gaunt:2012dd a factor $1/2$ is missing on the r.h.s.

[TABLE]

Just as graph 10b, graph c contributes to the evolution of $G$ in (94), giving a factor $\alpha_{s}\log(\mu_{y}^{2}/\Lambda^{2})$ at leading order. There is however no contribution to the evolution of $F_{\text{int,pt}}$ from $\mu_{y}$ to $\mu_{2}$ from this graph, so that in the cross section we have

[TABLE]

which again reproduces the result (97) obtained in momentum space, including the overall factor.

A fully analogous discussion can be given for the graph in figure 5b. One finds the same type of DGLAP logarithms as in (97), which is again reproduced correctly in the DPS cross section formula. Note that in this case the additional gluon is a virtual correction rather than emitted into the final state – a somewhat unusual situation in the context of DGLAP logarithms.

In conclusion, higher-order graphs with strong ordering of transverse momenta produce DGLAP type logarithms in 2v1 graphs, which can depend on $Q_{\text{min}}/\Lambda$ or $Q_{2}/Q_{1}$. Our expression for the regularised DPS cross section correctly reproduces these logarithms if we make an adequate choice of DPDs and scales. Specifically, we must

• •

use DPDs that have the correct short-distance behaviour, as discussed in section 3.3. Fixed-order perturbative results for this behaviour can be used reliably at scales of order $\mu_{y}$, as in (90) and (94).

• •

take scales $\mu_{1},\mu_{2}$ and $\nu$ in $\sigma_{\text{DPS}}$ as specified in (82), with possible modifications by factors of order unity.

The situation is quite different for 1v1 graphs. In the DPS region

[TABLE]

the double box in figure 11a develops a squared DPS logarithm in $d\sigma/dq_{1}^{2}$,

[TABLE]

which disappears however when one integrates over $q_{1}$ up to $Q_{\text{min}}$. This is inevitable since the integrated 1v1 cross section depends on a single dimensionful quantity $Q_{\text{min}}$ and hence cannot have any large logarithm. An infrared cutoff $\Lambda$ in the $q_{1}$ integration does not play any role here since the cross section is dominated by large $q_{1}\sim Q_{\text{min}}$, where strong ordering as in (99) is not possible.

One finds that the higher-order graph in figure 11b does have a DGLAP logarithm in the strong-ordering regime

[TABLE]

due to the additional loop integral $\int d\ell^{2}/\ell^{2}$. But the region (101) gives only a power suppressed contribution to the cross section integrated over $q_{1}$, where DPS kinematics is not dominant.888This point was overlooked in the discussion of DGLAP logarithms in Diehl:2016khr ; Diehl:2016rqr .

There is an exception to the statement just made. If at least one of the momentum fractions $x_{1},x_{2}$ is very small and the corresponding scale $\mu_{1},\mu_{2}$ is sufficiently large, then the evolution of $F_{\text{spl,pt}}$ from $\mu_{y}$ to $\mu_{1}$, $\mu_{2}$ can significantly flatten the $1/y^{2}$ behaviour of the initial condition (90). The underlying physics is radiation of soft gluons, which strongly enhances parton distributions at low $x$, as is well known from the evolution of ordinary PDFs. The enhancement by evolution of $F_{\text{spl,pt}}$ increases with the distance between $\mu_{y}$ and $\mu_{1},\mu_{2}$, i.e. it increases with $y$. In the 1v1 (and also the 2v1) contribution to the DPS cross section, this enhances the importance of the region $y\gg 1/Q$, where the DPS approximation is valid. In this region, the evolution of the DPDs correctly describes the DGLAP logarithms generated by graphs like the one in figure 11b. In the numerical studies presented in section 9.2, we find that this happens for $F_{gg}$, $F_{qg}$ and $F_{\bar{q}g}$, as well as for $F_{q\bar{q}}$ when one momentum fraction is much smaller than the other.

6.2 Dependence on the cutoff scale

As already discussed in section 4.1, our formalism relies on the fact that in the region of small $y\sim 1/Q$ the subtraction terms in the cross section (45) cancel against the DPS contribution, up to corrections that are beyond the accuracy of the calculation. Such corrections include higher-order terms in $\alpha_{s}$, whose parametric size we now investigate.

Let $F_{\text{fo}}(x_{i},\boldsymbol{y};\mu)$ denote the short distance approximation of a DPD, computed in fixed-order perturbation theory. It depends on $\mu$ via the scales of $\alpha_{s}$ and of PDFs or twist-four distributions (and at higher orders also via powers of $\log(y\mu)$ in the splitting kernels). In the cross section, this fixed-order approximation is used in two different ways:

For the small-$y$ limit of the DPDs used in the DPS cross section, $\sigma_{\text{DPS}}$, one takes $F_{\text{fo}}(x_{i},\boldsymbol{y};\mu_{y})$ as initial condition and evolves it to the scales $\mu_{1}\sim Q_{1}$ and $\mu_{2}\sim Q_{2}$ of the two hard scatters. 2. 2.

In the subtraction terms, $\sigma_{\text{1v1,pt}}$ and $\sigma_{\text{2v1,pt}}$, one takes $F_{\text{fo}}(x_{i},\boldsymbol{y};\mu_{h})$ evaluated at the same scale $\mu_{h}$ that is used in the hard scattering terms, $\sigma_{\text{SPS}}$ or $\sigma_{\text{tw4}}$, so as to cancel the contribution from the region $y\gg 1/Q$ order by order in perturbation theory.

The scale $\mu_{h}$ may be a combination of $Q_{1}$ or $Q_{2}$ or an independent third scale. For now we assume that $\alpha_{s}\mskip 1.5muL\ll 1$, where

[TABLE]

The case where this is not satisfied because the hard scales are widely different is discussed in the following subsection.

Consider now a calculation at NkLO ($k=0$ for LO), using DPDs evolved with NkLO kernels $P$ and NkLO expressions for the DPS subprocess cross sections $\hat{\sigma}_{1,2}$ and the short-distance approximations $\smash{F_{\text{fo}}}$. We have

[TABLE]

and its analogue for $\hat{\sigma}_{2}$. Here $\otimes$ denotes the usual Mellin convolution for the indicated variable, and we omit labels and sums for parton species. To understand (103) we observe that the l.h.s. would be exactly zero if $\hat{\sigma}_{1}$ and $P$ were calculated to all orders; the convolution of $\hat{\sigma}_{1}$ with the relevant parton densities would then be $\mu_{1}$ independent. With NkLO truncations on the l.h.s., the first term beyond the accuracy of the calculation is of order $\alpha_{s}^{k+1}\mskip 1.5mu\hat{\sigma}_{1}$. Furthermore, the highest power of $\log(\mu_{1}/Q_{1})$ in $\hat{\sigma}_{1}$ is $k$ according to the general pattern for logarithms related with the renormalisation or factorisation scale. In a similar fashion we can derive

[TABLE]

using that $k$ is the highest power of $\log(y\mu)$ in the NkLO approximation, according to (27). Finally, it is easy to see that for a quantity satisfying

[TABLE]

one has

[TABLE]

where the integration of the differential equation (105) between $\mu_{a}$ and $\mu_{b}$ gives an extra logarithm on the r.h.s. We thus obtain for the small-$y$ contribution to the 1v1 cross section

[TABLE]

using our initial condition $F_{\text{spl,pt}}=F_{\text{spl,fo}}$ at scale $\mu_{y}$. The longitudinal momentum arguments and convolutions are as in (2). The first line of the last expression in (6.2) is just the double counting subtraction term. We thus find

[TABLE]

with $L$ defined in (102). In the second step we have used that the full SPS cross section is of the same order of magnitude as the DPS subtraction term, regarding its power behaviour in $Q^{2}$ as well as the relevant coupling constants and PDFs. In full analogy one can derive

[TABLE]

Notice that the difference $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ is by construction dominated by short distances $y\sim 1/Q$, whilst the individual terms receive contributions from a wide range of $y$. We thus find that, for both 1v1 and 2v1 contributions, the DPS term cancels against the double counting subtractions in the region $y\sim 1/Q$, up to corrections of relative order $\alpha_{s}^{k+1}$. This is beyond the accuracy of an NkLO calculation, as it should be. The only logarithms multiplying $\alpha_{s}$ in this context are logarithms of the ratios between different hard scales, which at this point we must require to be of moderate size.

We can now also quantify the dependence of the overall cross section $\sigma$ in (45) on the scale $\nu$. With two choices $\nu_{1}\sim\nu_{2}\sim Q$ the difference $\sigma(\nu_{1})-\sigma(\nu_{2})$ is just of the size given by (6.2) and (6.2). In analogy to the variation of renormalisation and factorisation scales, one may thus use the variation of $\nu$ to estimate uncalculated higher-order corrections.

An important consequence of (6.2) is that

[TABLE]

where we used that the subtraction term $\sigma_{\text{1v1,pt}}$ is dominated by $y\sim 1/Q$, given the $1/y^{2}$ behaviour of $F_{\text{spl,fo}}$. This means that we can use the $\nu$ variation of the 1v1 contribution to $\sigma_{\text{DPS}}$ to estimate the size of the subtraction term $\sigma_{\text{1v1,pt}}$ and thus of the SPS contribution to the cross section at the corresponding order in $\alpha_{s}$. This is of great practical value, given the considerable difficulty to compute SPS at this order. We will come back to this point in section 9.2.

6.3 Multiscale problems

The arguments leading to (6.2) and (6.2) are based on an expansion in $\alpha_{s}$ and break down when $\alpha_{s}\mskip 1.5muL$ is not small, for instance because there is a large hierarchy $Q_{1}\ll Q_{2}$ between the two hard scales. We now show that a modified choice of DPDs in the subtraction terms $\sigma_{\text{1v1,pt}}$ and $\sigma_{\text{2v1,pt}}$ can cope with such a situation.

Let us for definiteness describe the 1v1 subtraction; the 2v1 case can be treated in full analogy. We use the notation $F_{\text{spl,pt}}(x_{i},\boldsymbol{y};\mu_{b},\mu_{c}\mskip 1.5mu|\mskip 1.5mu\mu_{a})$ for the distribution that is initialised to

[TABLE]

and evolved according to the double DGLAP equations to the scale $\mu_{b}$ ($\mu_{c}$) for the parton with momentum fraction $x_{1}$ ($x_{2}$). The small-$y$ limit of the distributions used in $\sigma_{\text{DPS}}$ is $F_{\text{spl,pt}}(x_{i},\boldsymbol{y};\mu_{1},\mu_{2}\mskip 1.5mu|\mskip 1.5mu\mu_{y})$ and its analogue for $F_{\text{int,pt}}$. Let us now introduce a profile function $p(u;\mu_{a},\mu_{b})$ that interpolates monotonically between the limiting values $\mu_{a}$ at $u=0$ and $\mu_{b}$ at $u\to\infty$, satisfying

[TABLE]

We can then take distributions

[TABLE]

in $\sigma_{\text{1v1,pt}}$. These distributions tend to $F_{\text{spl,fo}}(x_{i},\boldsymbol{y};\mu_{h})$ for $y\gg 1/\nu$, so that the subtraction term cancels the contribution of $\sigma_{\text{SPS}}$ in this region. For $y\sim 1/\nu$ the distributions tend to $F_{\text{spl,pt}}(x_{i},\boldsymbol{y};\mu_{1},\mu_{2}\mskip 1.5mu|\mskip 1.5mu\nu)\approx F_{\text{spl,pt}}(x_{i},\boldsymbol{y};\mu_{1},\mu_{2}\mskip 1.5mu|\mskip 1.5mu\mu_{y})$, so that in that region the subtraction term cancels against $\sigma_{\text{DPS}}$. Distances $y\ll 1/\nu$ are of course removed in the cross section by the cutoff function $\Phi(y\nu)$. Using the same technique as in section 6.2, one can show that the cancellations just described are up to corrections of order $\alpha_{s}^{k+1}$ at NkLO, without any large logarithms. In particular, this ensures that the $\nu$ dependence of the overall cross section is beyond the accuracy of the computation.

As an alternative, one may replace $\nu$ in (113) with $\min(\mu_{1},\mu_{2})$, which requires less computation when $\nu$ is varied in the cross section (while $\mu_{1},\mu_{2},\mu_{h}$ are kept fixed). With this choice, the DPDs in the subtraction term reduce to $F_{\text{spl,fo}}(x_{i},\boldsymbol{y};\mu_{h})$ in the case of equal hard scales, $\mu_{1}=\mu_{2}=\mu_{h}$.

Profile functions for the renormalisation scale have been used in other contexts, see e.g. Ligeti:2008ac ; Abbate:2010xh . An suitable function for our purpose is

[TABLE]

where $u_{1}\sim 1$, $u_{3}\gg 1$ and $u_{2}={\tfrac{1}{2}}(u_{1}+u_{3})$. Functions that are piecewise polynomials have been used for instance in Abbate:2010xh ; Stewart:2013faa . To estimate theoretical uncertainties, one can vary the profile function $p$, for instance by varying $u_{1}$ and $u_{3}$. For a detailed discussion we refer to Stewart:2013faa .

6.4 Approximation of the intrinsic distribution

As discussed earlier, the correct small-$y$ limit $F_{\text{int,pt}}$ at leading order in $\alpha_{s}$ is obtained by making the ansatz (94) at scale $\mu_{y}$ and then evolving $F_{\text{int,pt}}$ to the scales $\mu_{1},\mu_{2}$ relevant for the DPS cross section. In practical terms, this requires an ansatz or model for the twist-four distribution $G$ at some reference scale $\mu_{0}$ and subsequent evolution of $G$ to the scale $\mu_{y}$. The technical implementation of twist-four evolution is rather involved, and one may wish to avoid it in situations where the corresponding loss of precision is tolerable.

In the model used later in this paper we indeed follow this path, taking as limiting form for small $y$

[TABLE]

at some fixed reference scale $\mu_{0}$ and then evolving $F_{\text{int,pt}}$ to the final scales $\mu_{1},\mu_{2}$. For the r.h.s. we will take a form proportional to the product $f(x_{1},\mu_{0})\,f(x_{2},\mu_{0})$ of ordinary PDFs (this model is restricted to DPDs in the colour singlet channel). Multiplying with a suitable function of $y$, one readily obtains a model for the intrinsic DPD at any value of $y$. This is indeed standard in the existing literature.

As follows from the discussions below (30) and below (94), evolution of the ansatz (115) to scales $\mu_{1},\mu_{2}$ correctly takes into account diagonal interactions as in figure 6c whilst neglecting non-diagonal interactions as in figure 6a and b, which in the correct treatment contribute to the evolution of $G$ from $\mu_{0}$ to $\mu_{y}$. For $\sigma_{\text{2v2}}$, these interactions are irrelevant since the small-$y$ region is power suppressed. Furthermore, we can expect them to be of limited importance for $\sigma_{\text{2v1}}$ in the situation described at the end of section 6.1. Evolution then leads to a $y$ dependence of $F_{\text{spl,pt}}$ that is significantly flatter than the $1/y^{2}$ behaviour obtained at fixed order. This shifts the dominant region of the $y$ integration in $\sigma_{\text{2v1}}$ to larger values, so that on the $F_{\text{int,pt}}$ side there remains little room for the evolution of $G$ between $\mu_{0}$ to $\mu_{y}$. The limited numerical importance of non-diagonal interactions for 2v1 diagrams in such situations was also observed in Gaunt:2012dd , where calculations were done in transverse momentum space.

7 Collinear DPDs in momentum space

As we argued in section 4, there are several advantages to work with DPDs in $\boldsymbol{y}$ space when setting up a factorisation formula for the overall cross section. We now consider collinear DPDs depending on the Fourier conjugate transverse momentum $\boldsymbol{\Delta}$, showing that they can be connected with our formalism in a useful way. The motivation for this is twofold.

• •

Sum rules for momentum space DPDs at $\Delta=0$ have been proposed in Gaunt:2009re , where it was stated that they are preserved by the leading-order inhomogeneous evolution equations in Kirschner:1979im ; Shelest:1982dg ; Snigirev:2003cq ; Gaunt:2009re ; Ceccopieri:2010kg . If it can be shown that they are valid at some starting scale, this will provide valuable (albeit nontrivial) constraints on DPDs. Of course, no such sum rules can hold for $\int d^{2}\boldsymbol{y}\,F(x_{i},\boldsymbol{y})$ since this integral is divergent at small $y$.

• •

Previous work on DPS, in particular on the perturbative splitting contributions Blok:2011bu ; Blok:2013bpa ; Ryskin:2011kk ; Ryskin:2012qx , is formulated in transverse momentum space. The results derived below will facilitate the comparison of that work with ours, which is the subject of section 8.

Both aspects only concern colour singlet DPDs, to which the rest of this section is restricted. Colour non-singlet DPDs have an additional rapidity dependence described by Collins-Soper equations, which are multiplicative in $\boldsymbol{y}$ space (see Buffing:2016wip ) and hence more complicated in momentum space.

The function $\Phi$ that regulates the DPS cross section in our formalism can also be used to introduce momentum space DPDs. We define

[TABLE]

omitting the superscript indicating the colour singlet channel for brevity. With the conditions on $\Phi$ stated in (• ‣ 4), the integral converges at $y=0$. Throughout this section, we restrict ourselves to the two choices of $\Phi$ in table 1 and assume that $\nu\gg\Lambda$. The dependence of $F_{\Phi}$ on $\nu$ is then given by

[TABLE]

The derivative of $\Phi(y\nu)$ strongly suppresses the nonperturbative region $y\sim 1/\Lambda$ in the integral (or removes it altogether if $\Phi$ is a step function). This allows us to replace $F(x_{i},\boldsymbol{y},\mu)$ with its perturbative splitting part on the r.h.s. The contribution from $F_{\text{int,pt}}(x_{i},\boldsymbol{y},\mu)$ is suppressed by $\Lambda^{2}/\nu^{2}$. Inserting the expressions (3.2) or (3.2) for $F_{\text{spl,pt}}(x_{i},\boldsymbol{y},\mu)$ at $\epsilon=0$ one obtains

[TABLE]

in general, where

[TABLE]

For brevity, we have omitted labels and sums for parton species and polarisations.

For $\Delta=0$ we have $I_{0}=1$ and recognise in (118) the inhomogeneous term of the evolution equation for a conventional momentum space DPD $F(x_{i},\boldsymbol{\Delta},\mu)$ defined by $\overline{\text{MS}}$ subtraction. At leading order, our momentum space DPD $F_{\Phi}(x_{i},\boldsymbol{0},\mu,\mu)$ thus obeys the same evolution equation as $F(x_{i},\boldsymbol{\Delta},\mu)$. At small $\Delta$ the inhomogeneous term is multiplied by $I_{0}(\Delta,\mu,\mu)=1+\mathcal{O}(\Delta^{2}/\mu^{2})$. The higher-order version (119) at $\Delta=0$ has the same structure as the inhomogeneous term in the evolution equation proposed in Ceccopieri:2010kg , although the evolution kernel may not be the same.

Let us now compute the momentum space analogue of $F_{\text{spl,pt}}(x_{i},\boldsymbol{y},\mu)$, both with the cutoff regularisation (116) and with $\overline{\text{MS}}$ renormalisation of the splitting singularity. This gives the correct result for the full momentum-space DPD in the limit $\Delta\gg\Lambda$, when the exponent $e^{i\boldsymbol{y}\boldsymbol{\Delta}}$ suppresses large $y\sim 1/\Lambda$ in the Fourier integral, so that the perturbative splitting approximation can be used.

Restricting ourselves to leading order, we start with (3.2) and use the Fourier integrals given in (167) and (A). For $\Delta\ll\nu$ we find

[TABLE]

with

[TABLE]

To compute the usual momentum space DPD, we take the Fourier transform (9) of $F_{\text{spl,pt}}(x_{i},\boldsymbol{y})$ in $D=4-2\epsilon$ dimensions, subtract the UV divergence and then set $\epsilon=0$. This gives

[TABLE]

where we have abbreviated $x=x_{1}+x_{2}$ and $v=x_{1}/(x_{1}+x_{2})$. Here $C_{\epsilon}$ denotes the $\overline{\text{MS}}$ subtraction term for the splitting singularity. Using (165) we get

[TABLE]

where $P^{\prime}(v,\epsilon)=\partial P(v,\epsilon)/\partial\epsilon$. The $\overline{\text{MS}}$ counterterm is equal to

[TABLE]

so that we obtain

[TABLE]

We can use the preceding results to compute the matching between the DPD defined with a $\Phi$ cutoff and the one obtained by $\overline{\text{MS}}$ renormalisation. Since the integral in (116) is finite, it can be smoothly continued to $D=4-2\epsilon$ dimensions. We then have

[TABLE]

The factor $[1-\Phi(y\nu)]$ strongly suppresses the region $y\gg 1/\nu$, so that on the r.h.s. we can use the fixed-order perturbative splitting contribution, evaluated in $4-2\epsilon$ dimensions. This gives

[TABLE]

for all $\Delta$, and we have a smooth limit for $\Delta\to 0$. For $\Delta\ll\nu$ we can use (121) for $\Phi(u)=\Theta(u-b_{0})$ with (127) and obtain

[TABLE]

At $\Delta=0$ we can in particular achieve equality of the two DPDs by choosing the scale

[TABLE]

With the splitting kernel $P_{g\to q\bar{q}}$ from (24) this gives for instance a scale varying between $\nu=\mu\mskip 1.5mue^{-1/2}\approx 0.6\mskip 1.5mu\mu$ and $\nu=\mu$. These values respectively pertain to $v=1/2$ (i.e. $x_{1}=x_{2}$) and to $v=0$ or $1$ (i.e. $x_{1}=0$ or $x_{2}=0$).

8 Comparison with other work

Three groups have previously studied the consequences of the perturbative splitting mechanism on DPS cross sections, namely Blok et al. Blok:2011bu ; Blok:2013bpa , Ryskin and Snigirev Ryskin:2011kk ; Ryskin:2012qx , and Manohar and Waalewijn Manohar:2012pe . We refer to them as BDFS, RS and MW, respectively. In this section we compare the results obtained in our formalism with theirs. We restrict ourselves to collinear factorisation. To begin with, we rearrange our master formula (45) as

[TABLE]

and recall that one should take $\nu\sim\min(Q_{1},Q_{2})$ for the cutoff parameter. Let us briefly characterise the different terms in (132).

In all approaches just mentioned, the SPS cross section $\sigma_{\text{SPS}}$ is defined in the standard way. It is obtained by multiplying two PDFs with a hard scattering cross section and summing over all parton channels.

The term $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ does not have any equivalent in the work of BDFS and of RS. As discussed in section 4.1.2, this term is lacking the large $\log(Q/\Lambda)$ that is contained in $\sigma_{\text{DPS}}$ and built up by parton splitting taking place over a wide region $1/Q\ll y\ll 1/\Lambda$ in the 2v1 graph 1b. The approaches of BDFS and RS can thus only be valid at leading logarithmic accuracy. This limitation is explicitly stated in the work of BDFS Blok:2011bu . The same limitation holds of course in our approach if $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ is not computed.

The term $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$ includes the 2v2 and 2v1 parts of DPS, whereas the 1v1 part has its small $y$ approximation subtracted. As discussed at the end of section 6.1 and illustrated in section 9.2, there are situations where strongly ordered parton emissions in 1v1 graphs enhance the region $y\gg 1/Q$ for the initial parton splitting. In our approach, this effect is included via the evolution of the DPDs in $\sigma_{\text{DPS}}$ (which thus provides a small $x$ enhancement for $\sigma_{\text{2v2}}$, $\sigma_{\text{2v1}}$ and $\sigma_{\text{1v1}}$).

If such an enhancement does not take place, then the 1v1 part of DPS has an integrand going like $1/y^{4}$ at small $y$, so that the integral is dominated by the region $y\sim 1/Q$. In $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$ the contribution from this region then cancels up to higher orders in $\alpha_{s}$, as shown in section 6.2. In turn, the contribution from $y\sim 1/\Lambda$ is suppressed by $\Lambda^{2}/Q^{2}$ compared with $\sigma_{\text{SPS}}$ and thus of the same order as other power corrections beyond the accuracy of the calculation. In this situation, the dominant terms in $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$ are the 2v2 and 2v1 parts of DPS, provided that they have a small $x$ enhancement compared with generic power corrections to SPS.

8.1 The approach of Manohar and Waalewijn

The master formula for the overall cross section in the MW approach is given in eq. (14) of Manohar:2012pe . Adapted to our notation, it reads

[TABLE]

where $\hat{c}_{3}$ and $\hat{c}_{4}$ are hard scattering coefficients and $[\mskip 1.5mu\cdots\mskip 1.5mu]$ indicates UV renormalisation of the enclosed quantity. This means that the product $F(\boldsymbol{y})F(\boldsymbol{y})$ of DPDs is integrated over $\boldsymbol{y}$ in $2-2\epsilon$ transverse dimensions, and the resulting poles in $\epsilon$ are subtracted as usual in dimensional regularisation (we assume $\overline{\text{MS}}$ renormalisation, although this is not explicitly stated in Manohar:2012pe ). The terms going with $\hat{c}_{3}$ involve the renormalised product of a twist-two distribution in one proton and a twist-four distribution in the other. In the following, the terms with $\hat{c}_{4}$ and $\hat{c}_{3}$ in (133) will be called DPS and twist four terms, respectively. Via the perturbative splitting mechanism in one of the two protons, the operator in the DPS term mixes under renormalisation with the operators in the twist-four terms. This is discussed at LO in Manohar:2012pe , but we see no obstacle to extending the analysis to higher orders.

In section 4.1.2 we have seen that the computation of the hard scattering coefficient for graph 7 has a divergence from the region where the $g\to q\bar{q}$ splittings become collinear. In our formalism the term $\sigma_{\text{2v1,pt}}$ in (132) subtracts the contribution from this region. In the MW approach, the general subtraction formalism recalled in section 4 will yield an $\overline{\text{MS}}$ subtraction term for the collinear divergence, resulting from the renormalisation procedure described in the previous paragraph. For the choice $\Phi(u)=\Theta(u-b_{0})$ we see in (43) that $\sigma_{\text{2v1,pt}}$ is proportional to the $\overline{\text{MS}}$ expression $1/\epsilon+\log(4\pi)-\gamma$ plus an extra term $\log(\mu^{2}/\nu^{2})+P^{\prime}(0)/P(0)$, where $P(\epsilon)$ is the relevant splitting function in $4-2\epsilon$ dimensions. This extra term corresponds to the difference between the momentum space DPD $F_{\Phi,\text{spl,pt}}$ in (121) and its $\overline{\text{MS}}$ counterpart in (127). Up to this difference, the term $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ in our master formula (132) is thus equivalent to the twist four term in (133). The mismatch between the two versions may be understood as a scheme difference.

It remains to compare the DPS term in the MW approach with our term $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$. As the scale evolution of \bigl{[}\mskip 1.5mu\int d^{2}\boldsymbol{y}\,F(\boldsymbol{y})F(\boldsymbol{y})\mskip 1.5mu\bigr{]} derived by MW reflects perturbative splitting in one of the two protons, their DPS term should contain the contributions from 2v2 and 2v1 graphs. Our approach also contains the contribution of 1v1 graphs, with a subtraction of their small $y$ approximation. As discussed above, this part is negligible if DPS evolution effects are moderate. Both approaches should then give a valid representation of the overall cross section, with the possibilities to go beyond the leading logarithmic approximation and to systematically include higher orders in $\alpha_{s}$. On the other hand, strongly ordered emissions can significantly enhance the the region $y\gg 1/Q$ in 1v1 graphs. In our formalism, such emissions are resummed to all orders by the evolution of the DPDs, but we see no counterpart of this in the MW formulation.

MW argue that under renormalisation \bigl{[}\mskip 1.5mu\int d^{2}\boldsymbol{y}\,F(\boldsymbol{y})F(\boldsymbol{y})\mskip 1.5mu\bigr{]} does not mix with the product of two PDFs by pointing out that the 1v1 graph in figure 1a gives a quadratically divergent scaleless integral, which is zero in dimensional regularisation. This argument hinges on using the perturbative approximation for the $g\to q\bar{q}$ splitting at all values of $\boldsymbol{y}$ (or equivalently of the Fourier conjugate momentum $\boldsymbol{\Delta}$). This includes the infrared region, where the perturbative approximation does not work. It would be interesting to investigate this point further, but this is beyond the scope of the present paper.

8.2 The approaches of Blok et al. and of Ryskin and Snigirev

Since BDFS and RS use the same $\sigma_{\text{SPS}}$ as we do, and since they have no counterpart to our $\sigma_{\text{tw4}}-\sigma_{\text{2v1,pt}}$ (which can be neglected at leading logarithmic accuracy), it remains to compare our term $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$ with their expressions for the DPS cross section. In this comparison, we can take the leading logarithmic approximation of our result.

The approaches of BDFS and of RS are both based on separating the intrinsic and splitting contributions to a DPD at all values of $\boldsymbol{\Delta}$. As discussed in section 4, we avoid such a separation when formulating the factorisation formula for the cross section. When modelling the DPDs – which is of course unavoidable for phenomenology – it is however convenient to use this separation. In $\boldsymbol{y}$ space we thus make the ansatz

[TABLE]

for all $\boldsymbol{y}$. We demand that $F_{\text{int}}$ and $F_{\text{spl}}$ separately follow the homogeneous evolution equation (2) and that for $y\ll 1/\Lambda$ they tend to the correct short-distance limits $F_{\text{int,pt}}$ and $F_{\text{spl,pt}}$. At small $y$ the form (134) is hence unambiguous, whereas for large $y$ it represents a model. For $F_{\text{int,pt}}$ we make the approximation discussed in section 6.4. Then $F_{\text{int}}$ has a finite limit for $y=0$ rather than a weak singularity induced by non-diagonal twist-four evolution. Notice that $F_{\text{int}}$ and $F_{\text{spl}}$ now denote two terms in a DPD model, whereas in previous sections they were associated with different Feynman graphs in the discussion of factorisation.

BDFS and RS only consider DPDs for unpolarised partons in the colour singlet sector, neglecting spin and colour correlations in both the intrinsic and the splitting parts. We will therefore do the same when comparing to our formalism. We note that spin correlations tend to be washed out by evolution to high scales Diehl:2014vaa and that colour correlations are suppressed by Sudakov logarithms Mekhfi:1988kj ; Manohar:2012jr . For sufficiently large hard scales, the importance of such correlations will therefore be reduced.

In terms of the momentum space DPDs introduced in section 7, the DPS cross section at leading order in $\alpha_{s}$ is given by

[TABLE]

times the appropriate subprocess cross sections and a combinatorial factor. As was done in the papers cited above, we allow for different factorisation scales $\mu_{1},\mu_{2}$. We note that with the behaviour (122) of the splitting DPD at large $\Delta$, the integral (135) converges in the ultraviolet for both choices of $\Phi$. For definiteness, we take $\Phi(u)=\Theta(u-b_{0})$ in the following. We now analyse the 2v2, 2v1 and 1v1 contributions in turn.

The ultraviolet cutoff can be neglected in the intrinsic contribution:

[TABLE]

This is because $F_{\text{int}}(x_{i},\boldsymbol{y})$ is of generic size $1/\Lambda^{2}$ for $y\sim 1/\Lambda$ and has a smooth behaviour for $y\ll 1/\Lambda$. A cutoff $y>b_{0}/\nu$ in the integral therefore only removes a contribution of order $\Lambda^{2}/\nu^{2}$. If we make the same ansatz for the intrinsic DPD as BDFS or RS, we will thus obtain the same result for $\sigma_{\text{2v2}}$ as they do, up to power corrections that are beyond the accuracy of all all approaches discussed here.

For a typical ansatz, $F_{\Phi,\text{int}}(x_{i},\boldsymbol{\Delta})$ strongly decreases when $\Delta\gg\Lambda$. The 2v1 contribution to the integral (135) is therefore dominated by the region $\Delta\sim\Lambda$. In this region, the integral

[TABLE]

receives leading contributions from both small and large $y$. However, a large logarithm $\log(\nu/\Delta)$ is only built up in the region $1/\nu\ll y\ll 1/\Delta$, where we can use the perturbative form $F_{\text{spl,pt}}(x_{i},\boldsymbol{y})\sim 1/y^{2}$. To leading logarithmic accuracy, we can therefore compute the splitting contribution without a nonperturbative model, even for $\Delta\sim 1/\Lambda$, which we now do. We want to improve on our result (121) by including evolution effects. To do so, we start from the first order perturbative expression (3.2) at scale $\mu_{y}=b_{0}/y$ and then evolve using the Green function $D(x,\mu;\mu_{y})$, which solves the single DGLAP equation in the scale $\mu$ with the initial condition $D(x,\mu_{y};\mu_{y})=\delta(1-x)$. For simplicity we omit labels for and sums over parton species. We then have

[TABLE]

with

[TABLE]

The function $E(x_{i};\mu_{i},\mu_{y})$ depends weakly on $y$ due to logarithmic corrections from evolution and the running of $\alpha_{s}$. To make this explicit, we expand the quantities depending on $\mu_{y}$ (the Green functions, $\alpha_{s}$ and the parton density) around $\nu$. Keeping only leading logarithms, we obtain

[TABLE]

We show in appendix A that with such a $y$ behaviour, the Bessel function $J_{0}(y\Delta)$ in (138) can be replaced by an upper cutoff $\Theta(b_{0}-y\Delta)$, up to relative corrections in $\Delta^{2}/\nu^{2}$ and up to subleading logarithms $\alpha_{s}^{n+1}(\nu)\log^{n-2}(\nu/\Delta)$. We then have

[TABLE]

for $\Delta\ll\nu$. With the scale choice (82) our result (141) corresponds to the expressions of the perturbative splitting DPD in equation (18) of Blok:2011bu and in equation (8) of Ryskin:2012qx . Note that at leading logarithmic accuracy, the lower integration limit on $k$ in Blok:2011bu is equivalent to $\Delta$ in the region $\Delta\sim\Lambda$ that dominates the 2v1 cross section.

To leading logarithmic accuracy, the results for $\sigma_{\text{2v1}}$ obtained by both BDFS Blok:2011bu ; Blok:2013bpa and RS Ryskin:2011kk ; Ryskin:2012qx are therefore consistent with our approach (provided of course we make a suitable ansatz for $F_{\text{int}}$). In equation (6) of Ryskin:2012qx an upper cutoff $\min(Q_{1},Q_{2})$ is put on the $\Delta$ integration, but this is of no concern since as discussed above, the 2v1 term is dominated by $\Delta\sim\Lambda$.

The situation is different for $\sigma_{\text{1v1}}$, which RS compute using the form (141) we just derived for $\Delta\ll\nu$. The integral over $\Delta$ in the cross section then requires an upper cutoff, which RS choose as $\min(Q_{1},Q_{2})$, as we do for our cutoff $\nu$ in $y$ space. The dominant integration region for the 1v1 cross section is $\Delta\sim\nu$ in their approach and in ours, given that up to logarithms the splitting DPD is of order unity for both $\Delta\sim\Lambda$ and $\Delta\sim\nu$. For $\Delta\sim\nu$ our distribution $F_{\Phi,\text{spl}}$ differs from (141) even if we neglect evolution effects, which readily follows from the result (167) for the $y$ integral in (138). Furthermore, $F_{\Phi,\text{spl}}$ strongly depends on the choice of $\Phi$ in that region. This just highlights that $\sigma_{\text{1v1}}$ has no physical meaning by itself. In our formalism, the contribution from $\Delta\sim\nu$ in $\sigma_{\text{1v1}}$ is cancelled by the subtraction term $\sigma_{\text{1v1,pt}}$ , as it should be since the DPS approximation is not valid in that region at all.

The DPS cross section in the approach of BDFS is given entirely by 2v2 and 2v1 terms. To leading logarithmic accuracy, these terms agree with the ones in our approach. The term $\sigma_{\text{DPS}}-\sigma_{\text{1v1,pt}}$ in our approach has an additional contribution $\sigma_{\text{1v1}}-\sigma_{\text{1v1,pt}}$, which as already mentioned may be neglected unless DPD evolution enhances the region $y\gg 1/Q$ (or equivalently $\Delta\ll Q$) in $\sigma_{\text{1v1}}$. This possible enhancement is not captured by the BDFS approach.

In the RS approach, the DPS cross section includes 2v2, 2v1 and 1v1 terms, the first two of which agree with ours in the leading logarithmic approximation. If DPD evolution effects are moderate, then the 1v1 term of RS is dominated by large $\Delta\sim Q$, strongly depends on the upper cutoff imposed on the $\Delta$ integral, and leads to a double counting problem with $\sigma_{\text{SPS}}$ that is not solved in the RS approach. In the sense we have discussed in section 6.2, one may at best take the cutoff dependence of the 1v1 cross section as an estimate for the size of $\sigma_{\text{SPS}}$ and thus as an indicator for how serious the double counting problem is. Only if DPD evolution enhances the region $\Delta\ll Q$ so strongly that it dominates the 1v1 term does the RS approach give a valid representation of the overall cross section. The cutoff dependence of the 1v1 cross section is then reduced, as was emphasised in Ryskin:2011kk .

9 Quantitative illustrations

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

Convenient quantities for our purpose are double parton luminosities

[TABLE]

which directly appear in the DPS cross section. For definiteness, we use a sharp cutoff $\Phi(u)=\Theta(u-b_{0})$ to remove the UV region. As a model for colour singlet DPDs we take the sum of an intrinsic and a splitting part, as specified in (134) and in the paragraph following that equation. The luminosity (142) then naturally splits into separate terms for 2v2, 2v1 and 1v1. Colour non-singlet channels will not be considered, except in subsection 9.3, where they can be computed without recourse to a model.

9.1 Simplified analytic estimates

We shall see that DPD evolution has important quantitative effects on double parton luminosities. Nevertheless, let us start by deriving simple analytic expressions where evolution is neglected, but where the power behaviour, logarithms and other important factors can be easily identified.

For the intrinsic contribution we make the well-known product ansatz

[TABLE]

with a Gaussian form for the $y$ dependence. For the splitting contribution we take

[TABLE]

which is the perturbative expression at $O(\alpha_{s})$ times an exponential to suppress the large $y$ region. Sums and labels for the parton type are omitted here. For simplicity we use the same exponential in $F_{\text{int}}$ and $F_{\text{spl}}$. Since we ignore evolution for the time being, no scale dependence is taken into account in (143) and (144). To simplify even further, let us take equal momentum fractions $x_{i}=\bar{x}_{i}=x$ in both protons.

The parton luminosity for the 2v2 contribution then is

[TABLE]

As mentioned earlier, the UV regulator $\Phi(y\nu)$ only gives rise to a power suppressed effect. The 2v1 luminosity reads

[TABLE]

This is enhanced over $\mathcal{L}_{\text{2v2}}$ by $\log(\nu^{2}/\Lambda^{2})$ but suppressed by $\alpha_{s}/(2\pi)$. Furthermore, there are different factors regarding the parton densities. We note that the size of the splitting function strongly depends on the parton channel, as is exemplified by the cases $P_{g\to q\bar{q}}(1/2)=1/4$, $P_{q\to qg}(1/2)=10/3$ and $P_{g\to gg}(1/2)=27/2$. For the 1v1 luminosity we finally obtain

[TABLE]

We see the power enhancement over the 2v2 and 2v1 terms, as well as the quadratic dependence on the cutoff $\nu$.

It is easy to repeat the above exercise for our alternative regulator function $\Phi(u)=1-e^{-u^{2}/4}$. The integrals can be performed writing $\Phi(u)=u^{2}\int_{0}^{1/4}d\tau\,e^{-\tau u^{2}}$. For the leading term of the expansion in $\Lambda^{2}/\nu^{2}$, one obtains the same result for $\mathcal{L}_{\text{2v2}}$ as in (9.1), whereas in (9.1) and (9.1) one should replace

[TABLE]

respectively. The leading logarithm $\log(\nu^{2}/\Lambda^{2})$ in $\mathcal{L}_{\text{2v1}}$ is independent of the choice for $\Phi(u)$, but not the constant term. The corresponding ambiguity is compensated by the DPS subtraction term, as is the change of the numerical prefactor in $\mathcal{L}_{\text{1v1}}$.

9.2 Collinear Parton Luminosities

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

9.2.1 Setup

We use scale evolution at leading order in $\alpha_{s}$ and take equal scales $\mu_{1}=\mu_{2}=\mu$ for the two hard processes throughout. In this illustrative study, we fix the number of active quark flavours to $n_{f}=3$. This is sufficient for exploring the broad qualitative features of the DPDs and luminosities. Implementation of heavy quark flavours, as required for a realistic phenomenological investigation, would require some dedicated work on flavour thresholds and on memory management in our numerical code, which we postpone to future studies.

As in section 9.1, we write each DPD as a sum of intrinsic and splitting pieces:

[TABLE]

We initialise $F_{\text{int}}$ at a low scale $\mu_{0}=1\operatorname{GeV}$ using a product ansatz similar to (143):

[TABLE]

For the single PDFs $f_{a}(x,\mu_{0})$, we take the MSTW2008LO distributions Martin:2009iq . In (9.2.1) we have multiplied the PDFs by a function of the $x_{i}$ that does not affect the DPD at small $x_{i}$, but smoothly cuts it off near the kinematic boundary $x_{1}+x_{2}=1$. The function we use is that given in equation (3.12) of Gaunt:2009re , 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 Diehl:2014vaa , taking the width $h$ to be $x$ independent (corresponding to $h(x_{1},x_{2})$ of Diehl:2014vaa evaluated at $x_{1}=x_{2}=10^{-3}$) and setting each $h$ with $q^{-}$ indices to be the same as the one with $q^{+}$. Thus we have

[TABLE]

with

[TABLE]

For the splitting piece of the DPD we generalise our ansatz in (144), choosing an initialisation scale that goes to $b_{0}/y$ at small $y$ but freezes to a constant value $b_{0}/y_{\text{max}}$ when $y$ exceeds a value $y_{\text{max}}$ that marks the transition between perturbative and nonperturbative behaviour. This ensures that the single PDF and $\alpha_{s}$ in the splitting expression are never evaluated at too low scales. A suitable choice of scale is

[TABLE]

Such a prescription is very similar to the $b^{*}$ prescription used in TMD phenomenology Collins:1981va ; Collins:1984kg . Here we take $y_{\text{max}}=0.5\operatorname{GeV}^{-1}$, which is one of the values considered in the recent TMD study Collins:2014jpa . Using the same parton dependent Gaussian damping as in (9.2.1), we have

[TABLE]

The coupling $\alpha_{s}(\mu_{y})$ is determined by $3$ flavour running, starting with the MSTW2008LO value $\alpha_{s}(\mu_{0})=0.68183$. The single PDFs are obtained by taking the MSTW2008LO distributions at $\mu_{0}$ and evolving them according to the DGLAP equations for $n_{f}=3$.

To obtain the splitting and intrinsic DPDs at scale $\mu$, as in (149), the input forms just discussed must be evolved, starting from $\mu_{0}$ for $F_{\text{int}}$ and from $\mu_{y}$ for $F_{\text{spl}}$, according to the homogeneous double DGLAP equations. For this we use a modified version of the code developed in Gaunt:2009re . The modified code works on a grid in the $x_{i}$, $\mu$ and $\mu_{y}$ directions (the grid of the original code is in $x_{i}$ and $\mu$ only). The grid points in the $x_{i}$ directions are evenly spaced in $\log(x_{i}/(1-x_{i}))$, whilst those in the $\mu$ and $\mu_{y}$ directions are evenly spaced in $\log\mu$ or $\log\mu_{y}$. The integrals appearing in the double DGLAP equations are computed from points in the $x_{i}$ grids using Newton-Cotes rules (for details see Appendix A of Gaunt:2009re ), and evolution from one point in the $\mu$ grid to the next is carried out using the Runge-Kutta method.

The grid is set up with $88$ points in each $x_{i}$ direction, spanning $5\times 10^{-5}<x_{i}<1$, with $60$ points in the $\mu$ direction, spanning $\mu_{0}<\mu<170\operatorname{GeV}$, and with $60$ points in the $\mu_{y}$ direction, spanning $b_{0}/y_{\text{max}}<\mu_{y}<340\operatorname{GeV}$. According to the studies made in Gaunt:2009re , 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.

9.2.2 Numerical results

We begin with a study where the scale $\mu$ is equal to $Q_{1}=Q_{2}=80\operatorname{GeV}$ (as in the production of a $W$ boson pair). Taking the collider energy to be $\sqrt{s}=14\operatorname{TeV}$, we set $x_{1}$ and $\bar{x}_{1}$ in (142) to correspond to central production of the first system and $x_{2}$ and $\bar{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

[TABLE]

In figure 9.2.2 we plot $\mathcal{L}_{a_{1}a_{2}b_{1}b_{2}}(Y)$ in the range $0\leq Y\leq 4$ for the parton combinations $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}=u\bar{u}\bar{u}u+\bar{u}uu\bar{u}$, $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}=gggg$ and $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}=u\bar{d}\bar{d}u+\bar{d}uu\bar{d}$. The first parton combination appears e.g. in $ZZ$ production, the second is important in four-jet production, and the last appears in $W^{+}W^{+}$. For ease of language, we will refer to these parton combinations as the $u\bar{u}$, $gg$ and $u\bar{d}$ channels, respectively. We split the overall luminosity into contributions from 1v1 ($F_{\text{spl}}\times F_{\text{spl}}$), from 2v1 ($F_{\text{spl}}\times F_{\text{int}}+F_{\text{int}}\times F_{\text{spl}}$) and from 2v2 ($F_{\text{int}}\times F_{\text{int}}$). We vary $\nu$ by a factor of $2$ up and down around a central value of $80\operatorname{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 $\nu=80\operatorname{GeV}$, whilst the band is generated from the envelope of the functions with $\nu=40\operatorname{GeV}$, $\nu=80\operatorname{GeV}$ and $\nu=160\operatorname{GeV}$. In the case of the 2v2 contribution, the $\nu$ scale variation is negligible (as expected from basic considerations, see (9.1)), so this appears as a dashed line in each plot.

Figure 12: Double parton scattering luminosities $\mathcal{L}_{a_{1}a_{2}b_{1}b_{2}}(Y)$ for the production of two systems with $Q_{1}=Q_{2}=80\operatorname{GeV}$ at $\sqrt{s}=14\operatorname{TeV}$, one with central rapidity and the other with rapidity $Y$. The parton combinations $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}$ are $u\bar{u}\bar{u}u+\bar{u}uu\bar{u}$ (a), $gggg$ (b) and $u\bar{d}\bar{d}u+\bar{d}uu\bar{d}$ (c).

For the 1v1 contribution of the $gg$ and $u\bar{u}$ channels, we also plot a band generated by varying the prediction with $\nu=80\operatorname{GeV}$ by a factor $4$ up and down. This corresponds to a strictly quadratic cutoff dependence of $\mathcal{L}$ , i.e. to the variation of $\nu$ by a factor $2$ in the naive formula (9.1), 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 $u\bar{d}$ channel: there is no LO splitting process giving $u\bar{d}$, so that the scale variation (and central value) from (9.1) is zero in this case.

We immediately notice in figures 9.2.2a and b that the 1v1 contribution is generally much larger than the 2v2 and 2v1 contributions, and that it has an enormous $\nu$ 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 $\nu$ dependence of DPS. We also notice that at central rapidities in the $u\bar{u}$ channel, the 1v1 $\nu$ variation band essentially fills the band corresponding to a quadratic $\nu$ dependence, indicating small evolution effects in this case. By contrast, the 1v1 $\nu$ 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.

In figure 9.2.2c, the 1v1 contribution is small compared to the 2v2 piece and has a small $\nu$ dependence. This is because, as already mentioned, there is no LO splitting directly giving $u\bar{d}$ (generation of a $u\bar{d}$ pair requires at least two steps, such as $u\to u+g\to u+d+\bar{d}\,$). In this case, there is less of a need to compute $\sigma_{\text{SPS}}$ and the subtraction term $\sigma_{\text{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 $\alpha_{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 Blok:2013bpa ; Gaunt:2014rua at similar scales and $x$ values. In the $u\bar{u}$ channel, the magnitude of 2v1 and 2v2 is broadly the same but their shape is somewhat different. In the $u\bar{d}$ 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 $\nu$ variation in the $gg$ and $u\bar{u}$ channels is reduced by adding the remaining fixed-order terms in the cross section, even without considering a specific final state or having to compute $\sigma_{\text{SPS}}$. This is because only $\sigma_{\text{DPS}}$ and the subtraction terms depend on $\nu$. Just as we did for $\sigma_{\text{DPS}}$, we can factor the two parton-level cross sections out of $\sigma_{\text{1v1,pt}}$ and introduce a “subtraction luminosity” $\mathcal{L}_{\text{1v1,pt}}$. This is defined as in (142) except that the DPDs are replaced by the fixed-order splitting expression evaluated at scale $\mu$,

[TABLE]

For a given parton combination $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}$, one can directly subtract $\mathcal{L}_{\text{1v1,pt}}$ from the DPS luminosity $\mathcal{L}_{\text{DPS}}$ and compare the $\nu$ dependence of the result to that of $\mathcal{L}_{\text{DPS}}$ alone. Up to multiplication by the appropriate parton-level cross sections, the $\nu$ dependence of $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ is equal to that of the full cross section. On the other hand, the overall magnitude of $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ (which can be positive or negative) has no direct significance, since one needs to add the SPS contribution to obtain the physical cross section.

We plot $\mathcal{L}_{\text{DPS}}$ and $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ for the $u\bar{u}$ and $gg$ channels in figure 13. The lines are the values for $\nu=80\operatorname{GeV}$, whilst the bands show the variation when $\nu$ is varied by a factor of $2$ up and down. Note that for $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$, the curve with $\nu=80\operatorname{GeV}$ is always at the very top of the scale variation band. We see that the $\nu$ scale variation of $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ is indeed reduced compared to $\mathcal{L}_{\text{DPS}}$, with the reduction being much stronger for $u\bar{u}$ than for $gg$. The latter is consistent with our previous observation that evolution effects are weaker in the $u\bar{u}$ channel at central rapidities: if evolution effects are weak, $F_{\text{spl}}$ and $F_{\text{spl,fo}}$ have a similar $y$ dependence, so that $(F_{\text{spl}}\times F_{\text{spl}}-F_{\text{spl,fo}}\times F_{\text{spl,fo}})$ is flat in $y$ space and $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ varies weakly with $\nu$.

In section 6.2 we argued on generic grounds that the $\nu$ variation of the 1v1 contribution to $\sigma_{\text{DPS}}$ should be of the same order as the subtraction term $\sigma_{\text{1v1,pt}}$. This is confirmed in figure 13, where the size of $\mathcal{L}_{\text{1v1,pt}}$ (evaluated at central $\nu$) can be read off from the distance between the central curve for $\mathcal{L}_{\text{DPS}}$ and the top of the band for $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$. This finding allows us to sharpen the argument we already made in the discussion of figure 9.2.2: not only does a large $\nu$ variation of the 1v1 contribution to $\sigma_{\text{DPS}}$ indicate the need to include SPS and the associated subtraction term in the cross section, but the size of the $\nu$ variation may even serve as a rough estimate for the size of these terms.

Figure 14: Splitting DPDs $F_{\text{spl}}$ as a function of $y$, with $x_{1}=x_{2}=5.7\times 10^{-3}$ and $\mu=80\operatorname{GeV}$. The parton combinations shown are $u\bar{u}$ (a), $gg$ (b) and $ug$ (c). The vertical grey lines correspond to $y=b_{0}/(160\operatorname{GeV})$, $y=b_{0}/(80\operatorname{GeV})$ and $y=b_{0}/(40\operatorname{GeV})$.

To better understand the behaviour we have seen in the 1v1 luminosities for the $u\bar{u}$ and $gg$ channels, we now take a closer look at how evolution affects the $y$ dependence of splitting DPDs. In figure 9.2.2 we plot $F_{\text{spl}}(\boldsymbol{y};\mu)$ against $y$, setting $x_{1}=x_{2}=5.7\times 10^{-3}$ so that the $u\bar{u}$ and $gg$ DPDs are the ones used for the point $Y=0$ in figures 9.2.2 and 13. We also plot the $ug$ DPD, which mixes with $u\bar{u}$ and $gg$ under evolution. For comparison we show $F_{\text{spl,fo}}(\boldsymbol{y};\mu)$, as well as the initial condition $F_{\text{spl,ini}}(\boldsymbol{y})=F_{\text{spl}}(\boldsymbol{y};\mu_{y})$ given in (9.2.1), from which $F_{\text{spl}}(\boldsymbol{y};\mu)$ is obtained by evolution. The distributions are plotted in the range $b_{0}/(160\operatorname{GeV})<y<b_{0}/(40\operatorname{GeV})$, i.e. in the range over which the lower integration limit in $y$ is varied when we vary $\nu$ between $40\operatorname{GeV}$ to $160\operatorname{GeV}$. In this region, the exponential damping factor in our DPD model is irrelevant, so that $F_{\text{spl,ini}}$ and $F_{\text{spl,fo}}$ only differ by the scales taken in $\alpha_{s}$ and $f_{a_{0}}(x_{1}+x_{2})$.

For the initial conditions $F_{\text{spl,ini}}$ we note that the $u\bar{u}$ and $ug$ distributions are of similar size; the former is initialised by a larger PDF ($f_{g}$ instead of $f_{u}$) but has a smaller splitting coefficient $P(1/2)$ as we noted before (9.1). By contrast, the $gg$ distribution is much bigger; here both the initialising PDF and the splitting coefficient are large.

By construction, all three curves in each plot are equal at $y=b_{0}/(80\operatorname{GeV})$, when $\mu_{y}=\mu$. In all plots, $F_{\text{spl,ini}}$ has a more shallow $y$ behaviour than $F_{\text{spl,fo}}$. This difference is mainly driven by $\alpha_{s}$, as the PDFs do not vary so much between $\mu_{y}$ and $\mu$ at momentum fraction $x_{1}+x_{2}=1.14\times 10^{-2}$. One observes that the difference between $F_{\text{spl}}$ and $F_{\text{spl,ini}}$ is more significant for $gg$ and $ug$ than for $u\bar{u}$, i.e. that DPD evolution has a much stronger effect on the former two channels. This is to be expected at small $x$, since the $1/v$ behaviour of the splitting kernels $P_{gg}(v)$ and $P_{gq}(v)$ at small $v$ favours the radiation of a gluon with much smaller momentum than its parent parton, whereas the kernels $P_{qq}(v)$ and $P_{qg}(v)$ giving a quark stay finite for $v\to 0$.

An interesting point to note is that, whilst for $gg$ and $ug$ the curves for $F_{\text{spl}}$ are more shallow than for $F_{\text{spl,ini}}$, in the case of $u\bar{u}$ the $F_{\text{spl}}$ curve is actually steeper. The latter is surprising since — based on the experience with single PDFs at small $x$ — one may expect that forward evolution for $y>b_{0}/(80\operatorname{GeV})$ would always increase a DPD, and backward evolution for $y<b_{0}/(80\operatorname{GeV})$ would always decrease it. This indeed happens in the $gg$ and $ug$ channels, whilst forward evolution results in a decrease of the $u\bar{u}$ DPD. The reason for this is that $F_{ug,\mskip 1.5mu\text{spl}}$ and $F_{g\bar{u},\mskip 1.5mu\text{spl}}$, which directly feed into $F_{u\bar{u},\mskip 1.5mu\text{spl}}$ during evolution, are comparatively small. In the case of PDFs, $f_{g}$ is much larger than $f_{u}$ and hence can drive its small-$x$ evolution although the splitting function $P_{qg}$ does not favour the radiation of low-momentum quarks. The evolution of the $gg$ and $ug$ DPDs is driven by the large distribution $F_{gg,\mskip 1.5mu\text{spl}}$ and enhanced by the $P_{gg}$ splitting function.

Note that in all three channels, $F_{\text{spl}}$ has a smaller slope in $y$ than $F_{\text{spl,fo}}$. This implies that the $y$ integrand for the computation of $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ is positive for $y>b_{0}/(80\operatorname{GeV})$ and negative for $y<b_{0}/(80\operatorname{GeV})$. Therefore, the $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ curves for $\nu=40\operatorname{GeV}$ and $\nu=160\operatorname{GeV}$ lie below that for $\nu=80\operatorname{GeV}$. The curve for the central $\nu$ value thus lies at the top of the $\nu$ variation band, as we already observed in figure 13.

For $u\bar{u}$, the evolution from $\mu_{y}$ to $\mu$ turns out to have much the same quantitative effect as adjusting the scale in the fixed order expression from $\mu_{y}$ to $\mu$, such that $F_{\text{spl}}$ and $F_{\text{spl,fo}}$ end up extremely close together. This explains why at central rapidities the $\nu$ variation of the 1v1 $u\bar{u}$ contribution in figure 9.2.2a coincides almost exactly with the naive prediction from (9.1), which assumes that $F_{\text{spl}}$ depends on $y$ like $y^{-2}$ (as does $F_{\text{spl,fo}}$). Furthermore, if $F_{\text{spl}}$ and $F_{\text{spl,fo}}$ are very close to each other, then $\mathcal{L}_{\text{DPS}}-\mathcal{L}_{\text{1v1,pt}}$ is much smaller than $\mathcal{L}_{\text{DPS}}$, which we indeed see in figure 13a.

In the $gg$ channel, the modification of the $y$ slope by evolution is significant. In the range of figure 9.2.2, 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 $\mathcal{L}_{\text{DPS}}$ on the scale $\nu$ should then decrease. This was already anticipated in Ryskin:2011kk ; Ryskin:2012qx , where studies in the double leading logarithm approximation were performed.

Figure 15: DPS luminosities in the $gg$ channel with all $x$ values in the DPDs set equal. The collider energy $\sqrt{s}$ is fixed to $14\operatorname{TeV}$ (a) or to $100\operatorname{TeV}$ (b). The scales of the DPDs, set to $\mu=x\sqrt{s}$, are given at the top of the plot.

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 $\sqrt{s}$, so that the scale $\mu=Q_{1}=Q_{2}=x\sqrt{s}$ also varies with $x$. We make one plot with $\sqrt{s}=14\operatorname{TeV}$ and $5\times 10^{-4}<x<0.2$ (corresponding to $7\operatorname{GeV}<\mu<2800\operatorname{GeV}$) and another one with higher collider energy $\sqrt{s}=100\operatorname{TeV}$ and $10^{-4}<x<0.02$ (corresponding to $10\operatorname{GeV}<\mu<2000\operatorname{GeV}$). In our numerical code this requires a new DPD grid with larger $\mu$ and $\mu_{y}$ ranges, which was generated using 60 points in the $\mu$ and $\mu_{y}$ directions (as in the original grid). The resulting 1v1 luminosity is shown in figure 9.2.2, with bands for the $\nu$ variation and its naive version as in figure 9.2.2. For comparison we also show the 2v1 and 2v2 luminosities.

We see that the $\nu$ variation does indeed become progressively smaller compared to the central value as $x$ (and $\mu$) 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 $\nu$ variation becomes larger: for $\sqrt{s}=14\operatorname{TeV}$ the actual and naive $\nu$ bands essentially coincide towards the right of the plot, where values of the order of $x\sim 0.1$ are reached. Towards the left of each plot, where $x$ and $\mu$ become small, the 2v2 and 2v1 contributions begin to dominate over 1v1. For given $\mu$, this effect is more pronounced for $\sqrt{s}=100\operatorname{TeV}$ than for $\sqrt{s}=14\operatorname{TeV}$, since the $x$ values in the former case are smaller. We note that in the $u\bar{u}$ channel (not shown here) the actual and naive $\nu$ bands remain very similar throughout the kinematics of figure 9.2.2.

For the small $x$ values where the $\nu$ variation is small in figure 9.2.2, the $gg$ splitting DPD is significantly shallower in $y$ than the fixed-order $y^{-2}$ form. For $\sqrt{s}=14\operatorname{TeV}$ we find that around $x=2\times 10^{-3}$ it reaches $y^{-1}$, implying a logarithmic dependence on $\nu$ in $\mathcal{L}_{DPS}$, whilst a behaviour like $y^{-0.5}$ is reached around $5\times 10^{-4}$. For $\sqrt{s}=100\operatorname{TeV}$ the $y$ behaviour is like $y^{-1}$ around $x=4\times 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 $\mathcal{L}_{\text{DPS}}$ comes from large distances $y\gg 1/\nu$, where the DPS approximations used to derive $\sigma_{\text{DPS}}$ are valid. Also, the 1v1 part of $\sigma_{\text{DPS}}$ becomes clearly larger than its $\nu$ variation. As we argued earlier, the $\nu$ 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 $\mu=10\operatorname{GeV}$ in figure 9.2.2a. Using our argument that $\sigma_{\text{1v1,pt}}$ is of the same magnitude as the order of $\sigma_{\text{SPS}}$ that contains the lowest-order DPS-type loop, one may in this situation justifiably make predictions that include the DPS piece but omit the order of SPS just specified, as well as the associated subtraction term. This is encouraging, for instance in the context of four-jet production, where the computation of the relevant SPS order (namely NNLO) is well beyond the current state of the art. Note that lower orders in SPS should be computed and included, if possible.

Notice that when the $y$ behaviour of $F_{\text{spl}}$ becomes flatter than $y^{-1}$, the dominant $y$ region in the 1v1 part of $\mathcal{L}_{\text{DPS}}$ shifts to values $y\sim 1/\Lambda$, where one cannot compute the splitting DPD in perturbation theory and must rely on a model. Likewise, the 2v1 part of $\mathcal{L}_{\text{DPS}}$ becomes increasingly sensitive to the region $y\sim 1/\Lambda$ as soon as $F_{\text{spl}}$ becomes flatter than the fixed-order form $y^{-2}$.

Another kinematic regime where the 1v1 contribution to $\mathcal{L}_{\text{DPS}}$ becomes large compared with its $\nu$ variation and compared with $\mathcal{L}_{\text{1v1,pt}}$ is when the two hard systems have a large separation in rapidity. This reduction of the $\nu$ variation can be seen as $Y$ approaches $4$ in figure 9.2.2, but only in the $u\bar{u}$ and not in the $gg$ channel. In both channels, the effect becomes more pronounced once the rapidity separation of the hard systems is increased beyond $4$. To illustrate this, we make plots similar to figure 9.2.2 but now with one hard system at rapidity $Y$ and the other at rapidity $-Y$ (rather than one at $Y$ and the other at [math]), 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 $\nu$ variation in the 1v1 contribution is strongly reduced towards the right hand side of the plots, becoming hardly visible in the $u\bar{u}$ channel. Also notable is the fact that for large $Y$, the 2v2 contribution in this channel exceeds 1v1, which strongly decreases between $Y=0$ and $2$.

This reduction in $\nu$ 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\to 2$ splitting 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, $\mu_{y}$ and $\mu$, is increased. Thus, evolution again flattens the DPD compared to the naive $y^{-2}$ expectation and reduces the $\nu$ dependence of the DPS luminosities.

The effect is particularly drastic in the $u\bar{u}$ channel, because the lowest-order splitting $g\to u\bar{u}$ is inefficient at generating a pair with very different $x$ values (as $P_{g\to q\bar{q}}(v)$ goes to a constant at small $v$). Therefore, the repeated splitting mechanism described in the previous paragraph is strongly preferred, even though its last step is penalised by the lack of small-$v$ enhancement in $P_{gq}(v)$. In figure 18 we plot the $u\bar{u}$ DPD in the $x_{1}\ll x_{2}$ configuration relevant for $Y=4$ in figure 17a. It has a similar $y$ dependence as the $u\bar{d}$ DPD in equal kinematics, with the curve for $u\bar{d}$ crossing zero exactly at $y=b_{0}/\mu$ (by construction), and the curve for $u\bar{u}$ crossing zero rather close to this point. In the former case, the lowest-order $1\to 2$ splitting process is forbidden, whereas in the latter it is numerically almost irrelevant. The situation for $x_{1}\gg x_{2}$ is the same.

To end this section, we study polarised distributions and luminosities. We limit ourselves to the splitting part $F_{\text{spl}}$ of the DPD and to the 1v1 contribution to $\mathcal{L}_{\text{DPS}}$. In fact, we have little guidance for modelling the intrinsic part $F_{\text{int}}$, where a product ansatz as in (9.2.1) makes no sense. We note that in Diehl:2014vaa , different ansätze were made for polarised DPDs, and the effects of evolution on these distributions were investigated. It was found that, in many cases, evolution to higher scales leads to a suppression of the polarised DPD with respect to its unpolarised counterpart.

We initialise the polarised $F_{\text{spl}}$ at $\mu_{y}$ using the expression in (9.2.1) with the unpolarised $1\to 2$ splitting functions replaced by their polarised counterparts (given in Appendix B of Diehl:2014vaa ). These distributions are then evolved to the scale $\mu$ using the appropriate polarised double DGLAP equations (the required polarised splitting functions are collected in Appendix A of Diehl:2013mla ). The settings we used for the evolution code are the same as specified in section 9.2.1.

We consider the same scenario as in figure 9.2.2, i.e. the production of two systems with $Q_{1}=Q_{2}=80\operatorname{GeV}$ at $\sqrt{s}=14\operatorname{TeV}$, one with rapidity [math] and the other with rapidity $Y$, and now include polarised 1v1 contributions. In figure 19a, we reproduce figure 9.2.2a but include 1v1 luminosities for the polarised parton combinations $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}=\Delta u\Delta\bar{u}\Delta\bar{u}\Delta u+\Delta\bar{u}\Delta u\Delta u\Delta\bar{u}$ and $\delta u\mskip 1.5mu\delta\bar{u}\mskip 1.5mu\delta\bar{u}\mskip 1.5mu\delta u+\delta\bar{u}\mskip 1.5mu\delta u\mskip 1.5mu\delta u\mskip 1.5mu\delta\bar{u}$. For brevity we refer to these combinations as $\Delta u\Delta\bar{u}$ and $\delta u\mskip 1.5mu\delta\bar{u}$ in the following, recalling that $\Delta q$ and $\delta 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 Kasemets:2012pr (combinations like $u\bar{u}\mskip 1.5mu\Delta\bar{u}\Delta u$ also appear but are not shown here). To avoid cluttering the plots, we omit the unpolarised 2v1 luminosity, reminding the reader that it is of the same magnitude as 2v2. In figure 19b we repeat the exercise for the pure gluon channel, including the 1v1 luminosities for the polarised combinations $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}=\Delta g\Delta g\Delta g\Delta g$ and $\delta g\mskip 1.5mu\delta g\mskip 1.5mu\delta g\mskip 1.5mu\delta g$, referred to as $\Delta g\Delta g$ and $\delta g\mskip 1.5mu\delta g$, where $\Delta g$ and $\delta g$ respectively denote longitudinal and linear gluon polarisation.

We see that the 1v1 luminosities for $\Delta u\Delta\bar{u}$ and $\delta u\mskip 1.5mu\delta\bar{u}$ essentially coincide with the one for $u\bar{u}$ at central rapidities, with $\nu$ variation bands of very similar size. These observations can be explained by the fact that at $Y=0$ we initialise all three DPDs using the same expression up to a minus sign, with $P_{g\to q\bar{q}}=|P_{g\to\Delta q\Delta\bar{q}}|=|P_{g\to\delta q\mskip 1.5mu\delta\bar{q}}|$ at $v=1/2$, and that evolution has a rather weak effect on all three DPDs for central rapidities. This close agreement extends out to higher values of $Y$ for the $u\bar{u}$ and $\Delta u\Delta\bar{u}$ curves, owing to the fact that $P_{g\to q\bar{q}}=|P_{g\to\Delta q\Delta\bar{q}}|$ for general $v$. By contrast, the $\delta u\mskip 1.5mu\delta\bar{u}$ curve clearly falls below the other two for larger $Y$, because $P_{g\to q\bar{q}}>|P_{g\to\delta q\mskip 1.5mu\delta\bar{q}}|$ for $v\neq 1/2$.

In figure 19b, the $\delta g\mskip 1.5mu\delta g$ luminosity lies well below the one for $\Delta g\Delta g$, which in turn lies somewhat below the one for $gg$. This is mostly driven by the differences in initialisation expressions for the DPDs: at $v=1/2$ (relevant at $Y=0$) the relevant splitting function satisfy for instance $P_{g\to gg}:P_{g\to\Delta g\Delta g}:P_{g\to\delta g\mskip 1.5mu\delta 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 $\Delta g\Delta g$, which in turn should be roughly $50$ times bigger than the one for $\delta g\mskip 1.5mu\delta g$. This expectation is quite close to the actual luminosity ratios, which are $\sim 4$ and $\sim 60$, respectively. The remaining difference is due to evolution effects, which increase $gg$ rather considerably, increase $\Delta g\Delta g$ to a lesser extent, and hardly change $\delta g\delta 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 $\sigma_{\text{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 $\phi$, which couple to quarks with a Yukawa term $c\mskip 1.5mu\phi\mskip 1.5mu\bar{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 $n_{f}$. We will not compare the subtraction term to the full SPS contribution, but rather to the $gg$ initiated part of $\sigma_{\text{SPS}}$, which is the piece that contains the lowest order DPS-type $gg\to\phi\phi$ box diagram. This piece is gauge invariant and can hence be meaningfully considered by itself.

In keeping with the notation of the paper, let us denote the mass of each produced boson $\phi$ by $Q$. Let $Y$ be the rapidity of the diboson system in the $pp$ centre-of-mass frame, $\hat{s}$ its squared invariant mass, and $\beta=\sqrt{1-4Q^{2}/\hat{s}}$ the velocity of one boson in the diboson centre-of-mass frame. The $gg\to\phi\phi$ part of the SPS cross section can be written as

[TABLE]

where q_{\text{max}}=\beta Q\big{/}\sqrt{1-\beta^{2}} and the $gg\to\phi\phi$ matrix element squared, $|\mathcal{A}_{gg\to\phi\phi}|^{2}$, is averaged over the spin and colour of the incoming gluons. The momentum fractions $x$ and $\bar{x}$ in the gluon distributions are obtained from the constraints $\hat{s}=x\bar{x}s$ and $Y={\tfrac{1}{2}}\log(x/\bar{x})$.

The matrix element squared for $gg\to\phi\phi$ can be obtained from the $gg\to HH$ matrix element squared given in Glover:1987nx by making the replacement

[TABLE]

In particular, $|\mathcal{A}_{gg\to\phi\phi}|^{2}$ in (9.3) is related to the terms “gauge1” and “gauge2” in Glover:1987nx according to

[TABLE]

The right hand side of (159) is given in Glover:1987nx for general quark mass $m_{q}$. We evaluate this expression numerically for very small $m_{q}$, using analytic $m_{q}\to 0$ approximations where this is necessary to avoid numerical instabilities.

Now we turn to the subtraction term. The two elementary $q\bar{q}\to\phi$ cross sections in this term contain the kinematic constraint $2\pi\delta(x_{i}\bar{x}_{i}\mskip 1.5mus-Q^{2})$, which fixes the momentum fractions $x_{i}$ and $\bar{x}_{i}$ entering the two hard subprocesses at given $x,\bar{x}$ and $\beta$ up to a two-fold ambiguity. The DPS subtraction term is given by

[TABLE]

where we sum over all possible colour, spin, and quark number interference/correlation possibilities. The index $R$ on $F$ denotes the colour channel ($R=1$ for colour singlet, $R=8$ for colour octet), and in the sum over $a_{1}a_{2}\mskip 1.5mub_{1}b_{2}$ we sum over both unpolarised (e.g. $q\bar{q}\bar{q}q$) and polarised (e.g. $\Delta q\Delta\bar{q}\Delta\bar{q}\Delta q$) combinations. By $\{F\to I\}$ we denote the same expression, but with the quark number diagonal distributions replaced by the quark number interference ones (see section 2 of Diehl:2011yj ).

The splitting kernels for the relevant spin combinations are Diehl:2011yj

[TABLE]

and ${}^{8\!}P=-\,{}^{1\!}P/\sqrt{N_{c}^{2}-1}$, where $j,j^{\prime}$ 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

[TABLE]

where $j,j^{\prime}$ are the indices for transverse quark polarisation and where we have used the spin projection operators given in equation (2.90) of Diehl:2011yj .

Inserting (161) and (162) into (9.3), we obtain

[TABLE]

Note that both (9.3) and (163) contain a product of gluon PDFs evaluated at the same $x$ values, $g(x)\,g(\bar{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

[TABLE]

It is straightforward to show that $\Sigma$ is a function of the variable $\beta$ only. We plot $\Sigma$ for SPS and the subtraction term in figure 20, where the curve for $\Sigma_{\text{1v1,pt}}$ corresponds to the central choice $\nu=Q$ of the cutoff scale.

We see that $\Sigma_{\text{1v1,pt}}(\beta)$ is indeed generally of the same order of magnitude as $\Sigma_{\text{SPS}}(\beta)$. The agreement is perhaps surprisingly good, given that $\Sigma_{\text{1v1,pt}}(\beta)$ 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 $\beta\to 0$ and $\beta\to 1$, which correspond to the high energy limit and the threshold limit, respectively. It is to be expected that the agreement is especially bad at these points, since in the subtraction term we effectively assume that the integration over the squared transverse momentum $q^{2}$ of the scalar particles goes from zero to values of order $Q^{2}$. For both $\beta\to 0$ and $\beta\to 1$ this assumption becomes a poor one: near threshold the phase space in $q^{2}$ shrinks to zero, whilst in the high energy limit, $q^{2}$ can go up to size $\hat{s}\gg Q^{2}$.

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 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 $\Phi(y\nu)$ 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 transverse-momentum dependent DPDs and of collinear DPDs in colour non-singlet channels Buffing:2016wip .

The problem of double counting between DPS and other contributions — notably between DPS and SPS — is solved by subtraction terms as specified in (45) and (4.2), which are obtained in a simple way from $\sigma_{\text{DPS}}$ by replacing the DPDs with their appropriate short-distance limits. This paves the way for using the scheme at higher orders in $\alpha_{s}$, with calculations being considerably simpler for the subtraction terms compared with the full hard scattering process at the corresponding order. With a suitable choice of starting conditions and scales, specified in section 6.1, the DPS part of the cross section correctly resums DGLAP logarithms that are not included in the fixed order twist-four contributions.

Our scheme is naturally formulated with position-space DPDs $F(x_{i},\boldsymbol{y})$, but it is possible to relate the Fourier transform of $F(x_{i},\boldsymbol{y})\,\Phi(y\nu)$ to DPDs $F(x_{i},\boldsymbol{\Delta})$ 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 (130), and is somewhat similar to the matching between PDFs defined in different schemes such as the $\overline{\text{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 Blok:2011bu ; Blok:2013bpa and in Ryskin:2011kk ; Ryskin:2012qx to leading logarithmic accuracy.

For collinear DPDs, one can make a model ansatz consisting of two terms which in the limit $y\ll 1/\Lambda$ 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 $\nu$ of the 1v1 term in DPS provides an order-of magnitude estimate for the SPS contribution $\sigma_{\text{SPS}}$ (at the appropriate perturbative order), as it involves the same PDFs, overall coupling constants and kinematic region (small $\boldsymbol{y}$, corresponding to large transverse momenta and virtualities of internal lines). An alternative estimate is provided by the double counting subtraction term $\sigma_{\text{1v1,pt}}$, which by construction is dominated by small $\boldsymbol{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 $\beta$ of a boson in the boson pair rest frame is close to [math] 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 $\nu$. This demonstrates that, when including DPS in a cross section calculation, one must in general include $\sigma_{\text{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 $\nu$ 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 $\nu$ variation of the 1v1 DPS luminosity is considerably smaller than the central value. As we argued above, one may then justifiably neglect the order of $\sigma_{\text{SPS}}$ containing the first DPS-like double box, as well as the associated subtraction term, compared with $\sigma_{\text{DPS}}$. One scenario of this kind is when the flavour indices in both DPDs are $u\bar{d}$ (the luminosity with this parton flavour combination appears in $W^{+}W^{+}$ production). The suppression of the DPS-like double box in $\sigma_{\text{SPS}}$ is in this case related to the fact that perturbative splitting in the $u\bar{d}$ DPD starts at order $\alpha_{s}^{2}$ rather than at order $\alpha_{s}$. Further scenarios in which the $\nu$ variation in the 1v1 contribution is reduced are when $\sqrt{s}$ becomes very large compared to the hard scales $Q_{i}$, 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 $Q_{i}\ll\sqrt{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.

Appendix A Fourier integrals

In this appendix we collect a number of results for Fourier transforms in $2$ or $2-2\epsilon$ dimensions.

Fourier transform of a fractional power.

To compute the Fourier transform of $(\boldsymbol{k}^{2})^{-\lambda}$ in $2-2\epsilon$ dimensions we write $(\boldsymbol{k}^{2})^{-\lambda}=\Gamma^{-1}(\lambda)\int_{0}^{\infty}d\alpha\,\alpha^{\lambda-1}e^{-\alpha\boldsymbol{k}^{2}}$. One can then perform the Gaussian integration over $\boldsymbol{k}$, and subsequently the integral over $\alpha$. The result is

[TABLE]

As a corollary we obtain

[TABLE]

using $\boldsymbol{k}^{j}\mskip 1.5mue^{i\boldsymbol{y}\boldsymbol{k}}=-i\partial/(\partial\boldsymbol{y}^{j})\,e^{i\boldsymbol{y}\boldsymbol{k}}$.

Integrals involving $\Phi$.

We now compute the integrals in (48), which are also needed for (121) and (122), If $\Phi$ is the step function $\Theta(u-b_{0})$ then

[TABLE]

where $b_{0}$ is defined in (33) and ${}_{2}F_{3}(1,1;2,2,2;z)$ is a generalised hypergeometric function, which can be expanded as $1+\mathcal{O}(z)$ for small arguments. For $\Phi(u)=1-e^{-u^{2}/4}$ we can proceed as follows:

[TABLE]

where $E_{1}$ is the exponential integral. For small arguments one has $E_{1}(r^{2})=\log(1/r^{2})-\gamma+\mathcal{O}(r^{2})$. The corresponding integrals with $J_{2}$ instead of $J_{0}$ are given by

[TABLE]

and

[TABLE]

where in the second case we have used integration by parts. Both (169) and (170) behave like $1+\mathcal{O}(r^{2})$ for small $r$.

Connection between the Fourier-Bessel transform and a cutoff.

Consider the integral (138) with one term of the series (140). We will show that for $\Delta\ll\nu$

[TABLE]

The following argument is similar to the derivation given in appendix B of Bacchetta:2008xw . The integral on the r.h.s. is readily performed and gives

[TABLE]

The expression on the l.h.s. of (A) can be rewritten as

[TABLE]

using integration by parts. Since $J_{1}(u)=\mathcal{O}(u)$ at small $u$, the integral on the r.h.s. can be extended down to $u=0$ with an accuracy of $\Delta^{2}/\nu^{2}$. Rewriting the logarithm using the binomial series, we obtain

[TABLE]

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 $d_{k}$ whose values are not important here. Comparison with (172) gives the desired result (A).

Acknowledgements.

We gratefully acknowledge discussions with V. Braun, K. Golec-Biernat, D. Boer, A. Schäfer 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ödinger 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 Binosi:2003yf ; Binosi:2008ig .