Two-loop splitting in double parton distributions
Markus Diehl, Jonathan R. Gaunt, Peter Ploessl, Andreas Schafer
TL;DR
This paper calculates the next-to-leading order two-loop splitting functions for double parton distributions, providing essential ingredients for precise double parton scattering predictions at NLO.
Contribution
It presents the first two-loop matching calculations of DPDs onto PDFs and derives NLO 1->2 splitting functions for all partonic channels.
Findings
Computed two-loop matching of DPDs onto PDFs.
Derived NLO 1->2 splitting functions for all channels.
Verified sum rules and kinematic limits for the splitting functions.
Abstract
Double parton distributions (DPDs) receive a short-distance contribution from a single parton splitting to yield the two observed partons. We investigate this mechanism at next-to-leading order (NLO) in perturbation theory. Technically, we compute the two-loop matching of both the position and momentum space DPDs onto ordinary PDFs. This also yields the 1 -> 2 splitting functions appearing in the evolution of momentum-space DPDs at NLO. We give results for the unpolarised, colour-singlet DPDs in all partonic channels. These quantities are required for calculations of double parton scattering at full NLO. We discuss various kinematic limits of our results, and we verify that the 1 -> 2 splitting functions are consistent with the number and momentum sum rules for DPDs.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
DESY 19-024, CERN-TH-2019-013
Two-loop splitting in double parton distributions
Markus Diehl1, Jonathan R. Gaunt2, Peter Plößl3, Andreas Schäfer3
1 Deutsches Elektronen-Synchroton DESY, 22603 Hamburg, Germany
2 CERN Theory Division, 1211 Geneva 23, Switzerland
3 Institut für Theoretische Physik, Universität Regensburg, 93040 Regensburg, Germany
[email protected], [email protected], [email protected], [email protected]
July 18, 2019
Abstract
**Double parton distributions (DPDs) receive a short-distance contribution from a single parton splitting to yield the two observed partons. We investigate this mechanism at next-to-leading order (NLO) in perturbation theory. Technically, we compute the two-loop matching of both the position and momentum space DPDs onto ordinary PDFs. This also yields the splitting functions appearing in the evolution of momentum-space DPDs at NLO. We give results for the unpolarised, colour-singlet DPDs in all partonic channels. These quantities are required for calculations of double parton scattering at full NLO. We discuss various kinematic limits of our results, and we verify that the splitting functions are consistent with the number and momentum sum rules for DPDs. **
Contents
1 Introduction
Consider the integrated cross section for the production of some final state with associated hard scale at a hadron collider. As is well known, the dominant mechanism for this is single parton scattering (SPS), where one parton from each proton collides to form the state of interest. If the final state can be divided into two subsets and with associated hard scales and , then another mechanism is two separate parton-parton collisions, one yielding and the other yielding . This is called double parton scattering, or DPS. The integrated cross section for this mechanism is suppressed compared to that for SPS by , where denotes a hadronic scale and . Despite this, there are several motivations to study this process, and DPS occupies a special place amongst the power suppressed corrections for the production of . First, DPS may compete with SPS when the latter is suppressed by small or multiple coupling constants, with a well-known example being same-sign pair production [1, 2, 3, 4, 5, 6]. When followed by leptonic decays, this is an important process, since same-sign dilepton production is a search channel for new physics (see e.g. [7, 8, 9]). Second, although DPS is a power correction at the level of the integrated cross section, in differential cross sections there can be regions in which DPS is competitive with SPS. In particular, for small transverse momenta of and , DPS and SPS contribute at the same power [10, 11, 12, 13]. For several processes one also finds that DPS is an important contribution at large rapidity separation of and (see e.g. [14, 15, 7]). Third, the DPS cross section scales roughly as the fourth power of a parton distribution, whilst SPS scales as the second power (and other corrections scale as the third power). Thus for a given final state, as we increase the collider energy and probe smaller parton momentum fractions , DPS grows faster than the other contributions, and hence is especially relevant at the LHC and future, higher energy colliders. Finally, DPS is an interesting process to study in its own right: it provides information on the correlation between partons in the proton and thus reveals an aspect of hadron structure that remains inaccessible in single parton distributions.
A variety of DPS processes have been investigated at the ISR [16], the CERN SPS [17], the Tevatron [18, 19, 20, 21, 22, 23, 24, 25, 26], and the LHC [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 4] (see e.g. figure 4 of [35] and figure 15 of [38] for overviews). These measurements cover final states with jets, electroweak gauge bosons, heavy quarkonia, and open charm. The DPS contribution to pair production of heavy gauge or Higgs bosons (such as or ) is estimated to be small, but may well become important in the high-luminosity phase of the LHC, given the overarching goal to probe electroweak symmetry breaking with highest possible precision. A detailed review of experimental, phenomenological, and theoretical aspects of DPS can be found in [39].
Either one or both parton pairs in a DPS process can arise from a single parton via a perturbative “” branching process. At next-to-leading order (NLO) and above, this branching may produce additional partons that are emitted into the final state, in addition to the two “active” partons that subsequently enter hard collisions. The process in which we have splittings in both protons overlaps with a loop correction to SPS, in the sense that a given Feynman graph corresponds to either DPS or SPS depending on the momenta of its internal lines. Likewise, the process where we have a splitting in just one proton overlaps with a “twist 2 twist 4” power-suppressed process. Developing a theoretical QCD framework for DPS that includes splitting but avoids double counting between SPS and DPS (and the twist 2 twist 4 mechanism) is quite nontrivial. First approaches to this problem are given in [13, 40, 41, 42]. More recently, some of us developed a new formalism [43] that possesses several advantages over the previous approaches. This involves a description of the DPS cross section in terms of so-called position space double parton distributions (DPDs) , where is the separation in transverse space between the two partons.
One of the advantages of the approach in [43] is that it can be formulated at all orders in , with higher-order contributions that can be practicably computed. This opens the perspective to performing DPS calculations at full next-to-leading order accuracy. The framework of [43] is designed to make maximal re-use of quantities that are already computed for SPS. The partonic scattering cross sections for DPS and for SPS are identical. NLO effects in parton emissions from the incoming legs are correctly described by including the known NLO corrections to the DGLAP splitting kernels [44, 45, 46, 47, 48, 49, 50, 51, 52] in the corresponding kernels of the “double DGLAP” equation, which describes the scale evolution of . A missing key ingredient, specific to DPS, is the effect of the splitting at NLO. This is the subject of the present paper. What is technically required in the notation of [43] is the two-loop kernel for the perturbative matching of the DPDs onto the PDFs at small .
Why should one consider performing such calculations? Apart from unknown higher-order corrections, a major source of uncertainty in DPS cross section calculations is the lack of knowledge of the nonperturbative “input” DPDs. Unlike ordinary parton distributions (PDFs), these are very poorly constrained due to the limited experimental information on DPS. One can expect knowledge of the DPDs to improve in future, as further information emerges from experiment, and as this information is combined with theoretical constraints [53, 54] and information emerging from lattice computations (see [55] for work in this direction). However, it is unlikely that the situation for DPDs will approach that for PDFs in the near-term future. Despite a considerable uncertainty on the shape of the nonperturbative DPDs in general, there are a number of motivations for analysing DPS at NLO accuracy. Concerning the partonic cross sections , it is well known from SPS that in many cases their LO expressions give only a rough approximation, especially for final states containing jets. On the side of DPDs, the splitting mechanism requires only PDFs as a nonperturbative input and can thus be brought under good theoretical control. Several studies have emphasised the quantitative impact of this mechanism on DPS [56, 57, 58, 59, 60, 61, 62, 43]. It is thus natural to go beyond the LO approximation for this contribution and compute the splitting at order . Moreover, for some parton combinations, this is the first order that gives a nonzero result. An example is the distribution, which is in particular relevant for same-sign pair production. The recent study in [63] finds that the splitting mechanism has some impact on observables even in such a situation. A general aim of going beyond LO is to get a handle on the perturbative convergence. As previously mentioned, favourable regions for DPS typically involve small values, and in this region there is potential that large logarithms can have a significant impact, and even spoil the perturbative convergence. By studying the effect of the NLO corrections on DPDs, and more generally on DPS luminosities and cross sections, we can get an indication of whether there are convergence issues caused (for example) by small logarithms. Finally, it is quite possible that by constructing appropriate ratios of cross sections (involving, for example, the same parton flavours but different scales), one may minimise the dependence on the nonperturbative input and more directly test the perturbative part of the cross section. In this case, NLO computations clearly should improve the overall accuracy.
In this paper we compute the matching coefficients at two loops for all partonic channels (i.e. for all combinations of parton labels ). We restrict ourselves to the matching coefficients for the colour-singlet, unpolarised DPDs. These are expected to be the dominant DPDs at small and high [64, 65]; in fact they are the only DPDs considered in most phenomenological analyses of DPS. The coefficients for the colour interference and polarised DPDs will be presented in future work.
Aside from DPDs in position space, one can also consider momentum space DPDs . The variable is Fourier conjugate to and represents the difference in transverse momentum between a parton entering the DPD operator in the amplitude and the corresponding parton in the conjugate amplitude. A description of DPS in terms of momentum space DPDs is considered in a number of works [12, 13, 56, 66, 61]. Furthermore, the momentum-space DPDs at are the objects that satisfy the number and momentum sum rules of [53], as was discussed in [54]. Similarly to the position space DPDs at small , momentum space DPDs at large can be matched onto PDFs with perturbative kernels . For the momentum space DPDs, the effect of splittings also shows up in the evolution equations, which contain an inhomogeneous term where the PDFs are convolved with “ evolution kernels” [67, 68, 66, 11]. In [54] it was shown that the kernels are related by sum rules to the usual DGLAP splitting functions at all orders in . Finally, one can introduce matching kernels that allow one to compute the distributions that correspond to a given set of position space DPDs. This allows one to make use of the sum rule constraints of [53] also if one calculates DPS cross sections in terms of position space DPDs. The computation of the kernels at the leading order level is given in section 7 of [43] and will be extended to NLO in the present work.
This paper is organised as follows. In section 2 we define the kernels , and , analyse their renormalisation scale behaviour, and derive relations between them. The general structure of the matching between momentum and position space DPDs is discussed in section 3. In section 4 we give details of the two-loop computation of the matching kernels. We present our results and discuss various properties of them in section 5. A brief summary of this paper is given in section 6.
2 General analysis: ultraviolet behaviour and scale dependence
In the next two subsections we briefly review the definitions of bare and of renormalised DPDs, setting up our notation on the way. In section 2.3 we then derive how renormalisation affects the splitting kernels for DPDs and their logarithmic dependence on the renormalisation scale . This will later allow us to obtain renormalised kernels from the computation of bare Feynman graphs for the relevant splitting processes.
2.1 Basic definitions and notation
To begin with let us recall that bare, i.e. unrenormalised, DPDs for partons and are defined by matrix elements
[TABLE]
where for quarks and antiquarks and for gluons. We use light-cone coordinates for any four-vector and write its transverse part in boldface, . It is understood that the transverse proton momentum is zero and that the proton polarisation is averaged over. The twist-two operators
[TABLE]
for unpolarised partons are constructed from bare field operators and are familiar from the definition of bare PDFs,
[TABLE]
In the present work, we consider colour-singlet DPDs, so that the colour indices of the quark or gluon fields in (2.1) are implicit and to be summed over. The form (2.1) holds in light-cone gauge , and in other gauges Wilson lines need to be inserted between the two parton fields. The DPDs in (2.1) depend on the transverse distance between the two partons. Distributions in transverse momentum representation are obtained by a Fourier transform
[TABLE]
where it is understood that bare quantities are evaluated in space-time dimensions. After ultraviolet renormalisation, which will be discussed in section 2.2, one obtains finite parton distributions and can set . Due to their different ultraviolet behaviour, renormalised DPDs in position and in momentum space are not Fourier conjugate to each other, unlike their bare counterparts.
In the limit of small or large , renormalised DPDs can be computed in terms of renormalised PDFs and perturbative splitting kernels as
[TABLE]
where is the renormalisation scale and we abbreviate
[TABLE]
with being the Euler-Mascheroni constant. The general form of the convolution integrals in (2.1) was already given in [43] and follows from boost invariance and dimensional analysis. The splitting mechanism expressed in (2.1) also induces an inhomogeneous term in the evolution equation for momentum space DPDs, which reads [54]
[TABLE]
where the homogeneous terms have the same form as in the evolution of and are given in (22) below. At leading order (LO) in , the kernel can be expressed in terms of the familiar LO DGLAP splitting functions, as specified in (37). The evolution equation (2.1) then takes the form that has been derived long ago [67, 68]. Furthermore, one has at LO, as shown in [11]. In the present work we compute all kernels , and at next-to-leading order (NLO) in .
We refer to and as “ splitting kernels” and to as “ evolution kernel”, where “” refers to having one initial parton that splits into the two partons specified in the DPD. Starting at NLO, the splitting process can involve further unobserved partons in the final state. At LO the flavour of is uniquely fixed for given , so that the sum over can be omitted in (2.1) and (2.1), but this no longer holds at NLO.
Let us now introduce some elements of notation that will be useful for the calculations in the following sections.
Dimensional regularisation.
Working in dimensions, we expand dependent quantities in a Laurent series around ,
[TABLE]
In particular, gives the finite piece and the residue of the pole.
Convolution integrals.
Following [54], we use a compact notation for the different types of convolution integrals that avoids giving explicit momentum fraction arguments. In the following, let be a function of whilst , , are functions of one momentum fraction only. We write
[TABLE]
where the integration region is determined by the support properties of the functions. If a one-variable function is a PDF then for and , and if a two-variable function is a DPD, then for , and . Corresponding support properties hold for renormalisation factors and splitting kernels (which may include delta functions and plus distributions at the endpoints of their support). In analogy to (9), we define and the combined convolution .
In (2.1) and in (2.1) we encounter a second type of convolution,
[TABLE]
Interchanging the order of integrations, it is straightforward to show that
[TABLE]
so that we can write these convolutions without the brackets. We also have
[TABLE]
where is the usual convolution for functions of a single momentum fraction. We will also write without brackets. One can rewrite (10) in the form
[TABLE]
which will be useful in section 5.
Inverse convolution (one variable).
We define as the solution of . We will only use this inversion for functions that have a series expansion in ; one then obtains a unique inverse order by order in perturbation theory. In particular, if , then .
Parton indices.
The notation for convolutions can be extended to include appropriate summation over indices denoting the parton type (quarks, antiquarks, gluons). For quantities and depending on a single momentum fraction, we write as usual
[TABLE]
Note that we do not use the summation convention for parton indices, i.e. sums are indicated explicitly when indices are given. For quantities and depending on two momentum fractions, we define
[TABLE]
For the combination of convolutions in the first and second momentum fraction, we then have
[TABLE]
2.2 Renormalisation and evolution of DPDs
Let us recall the renormalisation of DPDs, which is well-known at LO and has been formulated at arbitrary perturbative order in [54]. The twist-two operators in (2.1) contain short-distance singularities that are renormalised with a factor as111Our convention for renormalisation factors of composite operators corresponds to the one in [69, 70]. Other authors, such as the ones of [71], use instead.
[TABLE]
for ordinary PDFs. Correspondingly, is renormalised by a factor for each parton,
[TABLE]
In addition to the divergences of the twist-two operators, the momentum-space distributions have a splitting singularity, which arises from integrating over in (4) because at small . This singularity is renormalised additively with a factor as
[TABLE]
Defining the ordinary DGLAP evolution kernel by
[TABLE]
and using that bare matrix elements are independent, one obtains the usual DGLAP equations for by taking the derivative of (17). In analogy, one obtains the homogeneous double DGLAP equation
[TABLE]
from (18). Taking the derivative of (19), one obtains the inhomogeneous double DGLAP equation
[TABLE]
with the evolution kernel
[TABLE]
where was defined in the previous subsection. Making parton indices explicit, we have DPDs and PDFs . Correspondingly, the renormalisation factors and evolution kernels for PDFs read and as usual. For the quantities associated with splitting in DPDs we write and , omitting the subscript on and when parton indices are written out.
2.3 Splitting kernels at higher order
Let us now see how renormalisation affects the splitting kernels and . For unrenormalised DPDs at small or large , the splitting mechanism gives factorisation formulae
[TABLE]
with bare kernels and , where the prefactor of in is taken out for convenience.222Compared with [43], this prefactor includes an additional . According to (17) to (19), we then have
[TABLE]
with kernels
[TABLE]
for renormalised distributions.
implementation and coupling renormalisation.
For the bare momentum space kernel, we have a perturbative expansion
[TABLE]
where is the bare coupling. In the last step, we used the abbreviation (6) and introduced coefficients
[TABLE]
which are independent, because for dimensional reasons. As discussed in [54], we implement the prescription by dividing bare graphs of order by . The coupling renormalisation constant is then plus a sum of pure poles in . The standard choice in the literature is . An alternative choice was proposed by Collins [70], with the statement that for quantities for which ultraviolet divergences give at most one power of per loop, the two schemes give identical results for renormalised quantities at . We shall verify this explicitly up to order in section 2.3.1.
Up to order , we have
[TABLE]
where for brevity we do not display the dependence of and . Here , and is the number of active quarks. This gives
[TABLE]
for the running coupling in dimensions, whereas the renormalisation group derivative in dimensions reads
[TABLE]
Expanding (2.3) up to order , we obtain
[TABLE]
for the bare splitting kernel.
Renormalisation factors and splitting functions.
The perturbative expansion of the renormalisation functions and reads
[TABLE]
Note that has no tree level term. With the implementation of the scheme discussed below (2.3), all counterterms are pure poles in , and we have
[TABLE]
where the evolution kernels are expanded as
[TABLE]
The residues of in (34) are obtained by inserting (31) into (20) and (23) and then comparing the terms on both sides, taking into account that the evolution kernels are independent of , see e.g. [54]. Comparing the terms, one finds that and start at order . In particular, one has
[TABLE]
At LO, the inhomogeneous evolution kernel has the form
[TABLE]
with a kinematic constraint due to the fact that the splitting process gives exactly two final state partons at this order. The one-variable kernel is equal to , except that it has no terms and no plus prescription on factors. The convolutions in (36) hence turn into ordinary products:
[TABLE]
Note that for flavour diagonal transitions these expressions still contain distributions, which will be made explicit in (5) and (5).
2.3.1 Momentum space kernels
We expand the renormalised kernel (2.3) in dimensions as
[TABLE]
The coefficient equals the th order term of the bare kernel in (2.3), the coefficients with are products of lower order terms with renormalisation counterterms for the twist-two operators or the QCD coupling, whilst comes from the counterterm for the splitting singularity. Expanding in , we obtain the kernel in dimensions,
[TABLE]
with
[TABLE]
The absence of poles in yields the finiteness conditions
[TABLE]
where in all cases.
At LO we have
[TABLE]
where in the last step we used (34). has at most single poles, and the first condition in (42) implies
[TABLE]
The coefficients of the kernel (40) in dimensions thus read
[TABLE]
At NLO we have
[TABLE]
where in the second steps we used (34) again. The finiteness conditions (42) now read
[TABLE]
Using (36) and the fact the double poles in are obtained by replacing with , we get
[TABLE]
This relation fixes the double pole of in terms of the LO kernels and thus serves as a cross check of the two-loop calculation. For the single poles, we obtain
[TABLE]
from which we can extract the NLO evolution kernel . Collecting the finite contributions, we get
[TABLE]
where we used (45) to replace \bigl{[}W_{B}^{(1)}\bigr{]}_{0} by . Inserting (45) and (2.3.1) into the factorisation formula (25) for and taking the derivative, one obtains the inhomogeneous double DGLAP equation (22) up to order as a cross check.
To evaluate (2.3.1) one needs up to . The full expressions for in dimensions are easily computed along the lines of section 5 in [11] and read
[TABLE]
Here we introduced the factor
[TABLE]
with specified below (28). Since has the same kinematic constraint on and as in (37), the convolutions in (2.3.1) are products as in (2.3), with replaced as appropriate.
The square of the tree-level graphs for the splitting that give rise to (2.3.1) appears in many higher-order calculations and has been computed in many papers before. Note however that is computed for transverse momenta of and that differ by between the amplitude and its conjugate (see section 4.1). This is typically not the case in other contexts. Nevertheless, the expressions on the r.h.s. of (2.3.1) agree with the expressions for given in section 3 of [72].
Equivalence of scheme implementations.
Let us show that the two choices for in the implementation of the scheme specified below (28) give the same kernels in dimensions at LO and NLO. The choices coincide in the constant and the first-order term of the Taylor expansion of around , but they differ in the second-order term . Because has only a single pole in , the renormalised LO kernels are identical for the two choices, as are the associated renormalisation factors and . Since contains at most double poles in , the NLO renormalisation factor is the same for the two choices as well. The only remaining dependence on in the renormalised two loop kernel at can thus come from the terms
[TABLE]
where the dependence on follows from (28) and is due to the terms and in (2.3.1). Thanks to the last relation in (2.3.1), the dependence on cancels in this expression, which completes our argument.
2.3.2 Position space kernels
The bare kernels and are related by a Fourier transform according to (4). Using equation (E.1) in [73], one gets for the term in (2.3)
[TABLE]
with given in (6) and
[TABLE]
where denotes the Riemann function evaluated at integer argument . We define333Here and in (59) we use a convention that differs from equations (3.15) and (3.16) in [43].
[TABLE]
In contrast to (39) there is no term here, because there is no splitting singularity in space. By contrast, the counterterms for the twist-two operators and the QCD coupling are common to momentum and position space. One thus has
[TABLE]
For the kernel at we write
[TABLE]
The coefficients and obey exactly the same relations as and in (41) and (42), except that the sum over always starts at .
At LO and NLO, the relation (58) implies
[TABLE]
The kernel is finite and has at most single poles, so that one has only one finiteness condition,
[TABLE]
This is satisfied thanks to (60) for and , together with the third relation in (2.3.1). For the finite part of the kernel, the relations (60) imply
[TABLE]
where we used the relation (41) and its analogue for . This gives the following simple relation between the kernels in 4 dimensions:
[TABLE]
Explicitly, we have
[TABLE]
Inserting this into the factorisation formula (25) for and taking the derivative, one correctly obtains the homogeneous double DGLAP equation (21) up to order .
Higher orders.
One may wonder whether (2.3.2) and thus (63) holds at all orders . Consider therefore , where the term of contributes for the first time, namely to the finite part of but not to its poles. With (56) and (60) one has
[TABLE]
The sum multiplying is zero due to the constraint (42) for . We thus find that (2.3.2) remains valid for . One order higher, we have
[TABLE]
According to (42), the sum over m^{k}\,\bigl{[}W^{(4,m)}\bigr{]}_{-4} is zero for but not for . With (41) we get
[TABLE]
and thus find that (2.3.2) is no longer valid for .
3 Matching between momentum and position space DPDs
In the formalism developed in [43], DPS cross sections are calculated using position space DPDs. On the other hand, one of the few theoretical constraints we have on the form of DPDs is given by the sum rules introduced in [53], which are formulated for momentum space DPDs at . Due to the splitting singularity, position and momentum space DPDs are defined with a different ultraviolet renormalisation. As a consequence, they are not simply related by a Fourier transform but by a short-distance matching formula, which can be calculated in perturbation theory and was derived at LO in [43]. We now generalise this matching to higher orders in .
Consider a cut-off regularised momentum space DPD as defined in [43]:
[TABLE]
where is a function satisfying for and for with some . This function regulates the splitting divergence at small that one would encounter when naively Fourier transforming the physical position space DPD in dimensions. Let us first consider the case where is defined with a hard cut-off . The difference between the renormalised momentum space DPD (19) and (68) is then given by
[TABLE]
Here we have explicitly written out the splitting counterterm from (19) and indicated all dependence, using the shorthand notation in the last term. For brevity, the dependence of distributions and factors on is not shown. The relation (3) allows us to separate the matching into a part at that involves DPDs in dimensions and a dependent part that involves only DPDs at the physical point .
In following is understood that the cut-off scale is large enough to justify replacing with the perturbative splitting contribution (25), up to corrections in powers of . As discussed in section 3.3 of [43], the corrections of order arise from twist-three distributions, which are expected to be small at low and . The dominant power corrections are then of order .
We now turn to the case where is not a hard cut-off. The difference between two momentum space DPDs defined with different regulator functions and is given by
[TABLE]
This involves only at the physical point, because each regulator ensures that integrals over are finite in the ultraviolet. For distances of hadronic size, the expression in square brackets vanishes because in that limit, so that one can replace by its perturbative splitting form (25). Using (59) one can thus express the difference (70) in terms of and the integrals
[TABLE]
where the angular integration has been performed and has given rise to the Bessel function . It is thus straightforward to perform a conversion between DPDs defined with different regulator functions . In the following we will therefore limit ourselves to the choice .
3.1 Matching at zero
Let us first discuss matching at , which we wish to write as
[TABLE]
To compute the matching kernel , we replace in the first term on the r.h.s. of (3) by its perturbative splitting form (25). With the dependence in (57), this gives rise to the integrals
[TABLE]
Expanding as
[TABLE]
in analogy to (57), we thus find that the coefficients , , and are related as
[TABLE]
with from (56), where in the last step we have used the relation (58) between and . The coefficient originates from the splitting counterterm in (3) and thus reads
[TABLE]
In physical dimensions, i.e. at , the matching kernel can again be written as
[TABLE]
where and obey the same relations as and in (41). The coefficients fulfil finiteness relations analogous to those for in (42).
Since , the pole terms of and coincide at LO and NLO, which ensures the validity of the finiteness conditions for and fixes all coefficients of the logarithms in (77). For the nonlogarithmic term we obtain
[TABLE]
where in the second step we used (75), (76) and (56). The sum multiplying vanishes due to the constraint (42) for . In total, we thus have
[TABLE]
Higher orders.
Let us see whether the relation (79) remains valid at higher orders. To this end, we compute the coefficient
[TABLE]
where we have used (75) and (42). We thus find that (79) is no longer valid for . We can, however, derive from (75) a relation between the matching kernels and the position space kernels that is valid to all orders, namely
[TABLE]
In analogy to (2.3.2), we thus get
[TABLE]
which for is of course consistent with (2.3.2) and (79).
3.2 Matching at nonzero
For matching at nonzero , we need both terms on the r.h.s. of (3). Using the perturbative splitting form of (25) and the expansion (59) of the kernel , one readily obtains
[TABLE]
with corrections of and coefficients given by the integrals
[TABLE]
which arise after performing the angular part of the integration. We make the dependence explicit by writing this as
[TABLE]
where
[TABLE]
with . In the limits and , the above equations can be further simplified. In the former case the limiting form of (85) is easily obtained by Taylor expanding the Bessel function, yielding
[TABLE]
with numerical coefficients . The limiting behaviour for is obtained by writing out the polynomial series of \log^{\mskip 1.5muk+1}(z^{2}/a^{2})=\bigl{[}\mskip 1.5mu\log(\nu^{2}/\Delta^{2})+\log(z^{2}/b_{0}^{2})\mskip 1.5mu\bigr{]}^{\mskip 1.5muk+1} and extending the integration to infinity in (86). This gives
[TABLE]
with numerical coefficients . Plugging this back into (85) and rearranging the binomial sum, we obtain the following expression for the large behaviour:
[TABLE]
where the ellipsis denotes terms with fewer powers of logarithms.
To use the matching equation (83) truncated in , one must avoid large logarithms at higher orders. For the term with this requires one to take . Satisfying this condition, one obtains large logarithms in (89), so that the choice has to be avoided. On the other hand, no large logarithms appear for according to (87). Fixed order matching is thus possible both for and for . The latter choice is in fact unavoidable for matching at the point .
For a closed form of can be given in terms of the generalised hypergeometric functions , namely
[TABLE]
3.3 Scale independence of matching
We now discuss the dependence of the matching equation (83) on the cut-off scale and its behaviour under renormalisation group evolution. Let us first verify explicitly that the dependence is the same on the left and right hand sides. The logarithmic derivative of the l.h.s. is
[TABLE]
where in the last step we used the perturbative splitting form of in (25) and (59), as is appropriate for . For the r.h.s. of (83) one gets
[TABLE]
using the integral representation (84) of . With the relation (82) between and , we see that both sides of the matching equation have the same dependence at each fixed order in , as it should be. This statement holds up to power corrections in , due to replacing with its perturbative splitting form.
Let us now consider the renormalisation scale dependence of (83). Using the appropriate – homogeneous or inhomogeneous – double DGLAP equations, one obtains for the l.h.s.
[TABLE]
where the term with arises from the momentum space DPD. On the r.h.s. of (3.3) one can replace the difference of DPDs with the matching expression (83). Comparing this to the renormalisation group derivative of the r.h.s. of (83), one finds that the matching equation holds at all if
- •
obeys the inhomogeneous double DGLAP equation. At LO and NLO, this can be explicitly verified using the evolution equation for and the equality of and stated in (79),
- •
the dependent term
[TABLE]
satisfies the homogeneous double DGLAP equation. This follows from the fact that the splitting form
[TABLE]
of satisfies the homogeneous equation and that and satisfy the same differential equation in .
As is the case for the perturbative splitting forms of and , the double DGLAP equations are fulfilled only up to order if the perturbative series for the splitting or matching kernel is truncated at order .
4 Two-loop calculation
The aim of this work is to obtain the kernels , , , and introduced in sections 2 and 3 at NLO accuracy. To this end it suffices to calculate only the momentum space kernel , since from this , , and can readily be obtained as shown in the sections just mentioned. In this section, we describe the computation of the NLO momentum space kernel in some detail.
To begin with, we apply the factorisation formula (24) to the bare DPD of partons and in a parton . Expanding the formula in , we obtain for the second order term
[TABLE]
where denotes the bare PDF of parton in parton . In the second step we used the tree-level expression and the fact that because the corresponding loop integrals do not depend on any dimensionful scale. We thus obtain the bare two-loop kernel directly from the two-loop graphs for the bare DPD of partons and in an on-shell parton . We compute with massless quarks and gluons, so that the relevant graphs have both infrared and ultraviolet divergences, both being treated by dimensional regularisation.
4.1 Channels and graphs
Let us analyse at which order in a kernel is first nonzero. The conservation of quark and antiquark flavour implies that the splitting starts at
LO for , and , 2. 2.
NLO for , and , , ,
where and label quark flavours and may be equal or different. We will refer to these as “LO channels” and “NLO channels”, respectively. Further channels at the same orders are obtained by interchanging partons and or by charge conjugation, where the kernels for charge conjugated channels are identical. All other splitting processes start either at NNLO or N3LO.
Real emission graphs (or “real graphs” for short) are shown in figures 1 and 2 for LO and NLO channels, respectively. Further graphs for the same channels are obtained by complex conjugation and by interchanging the lines for partons and if . The unobserved parton radiated into the final state is uniquely determined for given , , and : it is a gluon for the LO channels and a quark or antiquark for the NLO channels.
For the LO channels, we have also virtual loop graphs (or “virtual graphs” for short). They are obtained by dressing the LO splitting graphs with one vertex correction or with one propagator correction for the lines of parton or (propagator corrections for the massless on-shell parton are zero in dimensional regularisation). This is depicted in figure 3.
At this point a brief comment is in order about the flavour structure for channels where and are only quarks or antiquarks. In graphs 2(k), 2(l) and 2(n) there are separate fermion lines for quarks and , which may or may not have the same flavour. We denote the corresponding kernels as specified in the figure. By contrast, graphs 2(m) and 2(o) are of a valence type in the sense that they involve a single quark flavour; we denote the corresponding kernels with a superscript . The complete kernels for a specific flavour transition are thus obtained as
[TABLE]
where .
Feynman vs. light-cone gauge.
We performed two independent calculations of the graphs making up the bare partonic DPDs – one in covariant Feynman gauge, and the other in light-cone gauge. For cut graphs in Feynman gauge, we use the Feynman rules given in appendix D of [73]. Not given there are the expressions for the four-gluon vertex. The expression to be used to the left of the cut can readily be found in textbooks (see e.g. section 3.2.4 in [70]), whereas to the right of the cut one needs to take the complex conjugate expression (just as for the quark-gluon vertex).
In addition to the graphs shown in figures 1, 2, and 3, one gets in Feynman gauge also graphs with Wilson lines attaching to the active partons, as already mentioned below (3). These Wilson line graphs can be obtained from the graphs without Wilson lines following the prescription illustrated in figure 4. Generally speaking, Wilson lines appear in two cases. The first is when an active parton, or , is a gluon emerging from a three-gluon vertex on one of the upper legs, as for instance in figures 1(a), 1(b), and 3(a). Such graphs give rise to two kinds of Wilson line graphs as shown in figure 4(a) and 4(b). In the second case, the active parton is a quark or antiquark originating from a quark-gluon vertex on one of the upper legs, as is e.g. the case in figures 1(h), 1(l), 3(f), and 3(g). This leads to only one Wilson line graph with a gluon attaching to a quark Wilson line, visualised in figure 4(c) and 4(d).
Finally, in Feynman gauge one has to include the corresponding Fadeev-Popov ghost version for all graphs containing closed gluon loops, such as in figure 3(a) and 3(d).
Kinematics.
For real graphs, both in LO and NLO channels, the general kinematic structure is illustrated in figure 5. We work in a frame where the incoming parton has plus momentum and zero transverse and minus momentum. The plus momenta and of the active partons and are equal in the amplitude and its conjugate, whereas their transverse momenta differ by as shown in figure 5. Minus momenta are assigned in the same way as the transverse ones. The parton that goes into the final state has plus momentum , and its momentum is uniquely fixed by the momenta of the active partons and , i.e.
[TABLE]
To obtain collinear momentum space DPDs, we have to integrate the expressions corresponding to the real graphs in figures 1 and 2 over , , and , which corresponds to having in the matrix element (2.1). As we are interested in collinear DPDs, we also have to integrate over the transverse momenta and , which is tantamount to the condition in (2.1).
The kinematic structure of the virtual graphs we encounter in the LO channels is exemplified in figure 6. Due to momentum conservation we now find that and , so that the integrations over and become redundant. However, there now is an additional loop momentum over which we have to integrate.
4.2 Performing the calculation
We now move on to outline the details of our calculation in Feynman and light-cone gauge. The initial computation of the diagram expressions is done using FORM [74] or FeynCalc [75].444In practice, FORM was used in the light-cone gauge calculation, whilst FeynCalc was used for the Feynman gauge calculation. Of course, either code could have been used for both computations. One then needs to perform the integrations over phase space and loop momenta. As mentioned before, there are two types of graph contributing to the bare partonic DPD at NLO: real emission graphs in which one parton is emitted into the final state, and virtual loop graphs with a loop in the amplitude or its conjugate and no emitted particles in the final state. We discuss the methodology for performing the integrations for the two cases in turn.
4.2.1 Real emission diagrams
For the real emission diagrams we have four different topologies. These are exemplified by figure 1(a), 1(b), 1(c), 1(d), and we refer to them as LD (lower diagonal), UD (upper diagonal), UND (upper non-diagonal) and T2B (top to bottom). The four-gluon graphs in figure 1(e), 1(f), and 1(g) can be considered to fit into the T2B topology, because all the denominator factors in these diagrams are covered by those in the T2B topology.
For each graph, we begin by performing the integration over minus components of the loop and phase space momenta. This is most easily done by changing variables as follows:
[TABLE]
The integral over is performed using the on-shell delta function, whilst the integrals over and can be performed using complex contour integration. For all topologies, one is able to close the integration contour on only one pole, both in the and in the integration. We note in passing that the variable transformation (99) is also useful in the general proofs of [76, 54].
As all plus components are fixed by the operators in (2.1), it remains to integrate over and at fixed . Every graph may be written in a form where its denominator factors are contained within the following set (this may require some relabelling of integration momenta, and use of the fact that the integrals are invariant under ):
[TABLE]
Graphs of all topologies have denominators and . In addition LD graphs have , UD graphs have and , T2B graphs have and , and UND graphs have and .
We express all numerator factors in the graphs in terms of the and , and then use integration by parts identities [77, 78], as implemented in LiteRed [79], to express the results for all graphs in terms of a small set of master integrals. For graphs with the topology UND the master integrals have the general form
[TABLE]
Specifically, we need the integrals
[TABLE]
as well as integrals related to these by the simple transformation .
For the remaining topologies, the master integrals are of the form
[TABLE]
and specifically we need
[TABLE]
as well as the integrals on the first line of (102) and integrals related to these by .
We apply the method of differential equations [80, 81, 82, 83] to compute these master integrals. The integrals depend on the external variables , and . Since is the only dimensionful external quantity, each master integral only depends on via an overall prefactor of a power of that is easily determined by dimensional analysis. The differential equations in thus only provide trivial information and are not used. We compute the master integrals by solving the differential equations in . For a particular master integral, solving the equation determines the integral up to a “constant” of integration, which here is actually a function of . We determine this constant of integration by computing the leading behaviour of the master integrals in the vicinity of the line . As a cross-check, we then verify that the obtained solution satisfies the differential equation in .
The system of differential equations has an almost entirely triangular matrix structure. For the integrals the system schematically looks as follows:
[TABLE]
where black squares denote entries of the form while black diamonds denote entries of the form . Here is the operator transforming a given master integral to the corresponding master integral with . This enables us to solve most of the system using elementary methods for first-order differential equations and forward substitution, working downwards from and . The two-denominator integrals can be solved exactly in dimensions, whilst for the remainder we obtain a series in up to the required order. The only instance in which we have a coupled system of differential equations is for and . Here we transform the system to the canonical basis [84] using Fuchsia [85], which allows us to obtain a solution for these two integrals as a series in . The system of differential equations for the integrals is fully triangular, so that application of elementary methods plus forward substitution suffices in this case.
The leading behaviour of the integrals for is computed using the method of regions [86]. Let us introduce the scaling parameter , and say that . Then for all master integrals the following region gives a leading contribution in :
[TABLE]
For only, we identify a further leading region, namely
[TABLE]
For each region, we use the appropriate scaling and drop terms in the denominator that are subleading in (so that the result is homogeneous in ). Following this approximation, every master integral has a sufficiently simple form to be solved to all orders in by the method of Feynman parameters. Then one adds together the contributions from the leading regions to obtain the leading behaviour in the limit . For any master integral, the region gives a non-integer power of for , namely . This is because all denominators behave like , whilst the phase space contributes . For the region gives an integer power of , namely , because behaves like , whilst all other denominators and the phase space behave like .
For each master integral, we checked our analytic results against a numerical computation of the integral performed by FIESTA 2 [87] at 10 randomly chosen points. The results were found to agree within the precision of the numerical computation.
When computing the expressions for the graphs, careful consideration of the limit is needed, because some graphs are singular in this limit. Fortunately, we know the full (all-order in ) behaviour of the graphs for from the boundary condition computation for the master integrals. For graphs associated with the non-UND topologies, this singularity is always regulated by the dimensional regularisation parameter , since the non-UND master integrals only have the leading region for . We make a distributional expansion of the relevant factors as follows:
[TABLE]
where we denote the plus distributions as
[TABLE]
Useful relations for these distributions are given in [88]. Our final results turn out to involve only and not with .
For graphs associated with the UND topologies, specifically in the contribution from , one can have terms which just go as and are not regulated by . This corresponds to a rapidity divergence and requires an additional regulator. In covariant gauge, the rapidity regulator is typically inserted via some modification of the Wilson lines, whilst in light-cone gauge it corresponds to a modification of the gluon propagator. As noted in section 3.5 of [11] or in [64], DPDs in the colour singlet channel do not suffer from rapidity divergences.555This is analogous to what happens in the single scattering sector, where for the collinear PDFs we have no rapidity divergences [89]. Therefore, the divergences associated with must cancel when we sum over graphs, including the virtual loop graphs.
There are various possible choices for the rapidity regulator, such as the analytic regulator of [90], the exponential regulator of [91], the “pure rapidity regulator” of [92], the regulator of [93], the CMU regulator of [94, 95], or the use of Wilson lines tilted away from the light cone as described in [70]. We tried each of the latter three options, and we find that when consistently implementing the regulator in all real emission and virtual loop graphs, the sum over graphs yields the same rapidity-divergence-free result for all regulators. Note that we always remove the rapidity regulator before performing ultraviolet renormalisation and setting , as is prescribed when using such regulators. Note also that for each regulator we have a distributional expansion analogous to (108). We will give more detail in a forthcoming publication devoted to two-loop results for the colour non-singlet channels. In this case, the rapidity divergences only cancel when one combines the bare DPDs with the appropriate soft factor.
4.2.2 Virtual loop diagrams
Just as for the real graphs, the integrations over minus momenta (including ) in virtual graphs are easily done using complex contour integration. The integrations over the transverse momenta and were performed using different methods in the Feynman gauge and light-cone gauge calculations. In the former case we again used integration by parts reduction, as we did for the real emission graphs. One can even re-use the same master integrals, namely and , where the latter is obtained from by interchanging . In the calculation using light-cone gauge, the adopted strategy closely followed that of [52] – indeed we were able to directly re-use many of the results given in Appendix A of that work.
Finally, we have to integrate over a finite range. This finite range is a result of the earlier integrations over minus components: when is outside a certain range, all poles move into one half-plane for one of these integrations, and the result is zero. The integration may be performed using standard integration techniques. However, this again requires a consistent treatment of endpoint singularities. We have to perform a distributional expansion analogous to (108) in cases where these singularities are regulated by , and when this is not the case we have to insert a rapidity regulator. In light-cone gauge, the graphs containing such a rapidity divergence are vertex correction graphs in which the line with momentum is a gluon, such as figure 3(f). In Feynman gauge the rapidity divergence is contained in the corresponding Wilson line graphs. Here we insert the same regulator as used in the real emission calculations.
5 Results
In this section, we present our results for the kernels and up to NLO and discuss some of their properties. Results for the kernels and can readily be obtained from the relations (2.3.2) and (79).
Support and singularity structure
Let us briefly discuss the singularity structure of the NLO results. This differs between the LO and NLO channels defined at the beginning of section 4.1. The kernels for LO channels have the general form
[TABLE]
where is either or . The kernels for NLO channels contain only a term . The functions and are regular in the region defined by the conditions , and , and is regular in the region . Note that these regions exclude the points and , where some kernels have power-law divergences of the form with integers and . These points are not reached in any physical double parton scattering process, where both and must be strictly positive. We note that the decomposition (110) is invariant under the simultaneous replacement
[TABLE]
where is a regular function. This leaves a freedom of choice in the form of and for a given kernel.
Inserting the form (110) into the convolution with a PDF and using the representation (13), one obtains one-variable distributions and , so that
[TABLE]
with
[TABLE]
The term with in (5) has the same structure as the convolution of the LO kernels or with .
Before presenting the kernels in detail, we note that
[TABLE]
Analogous relations hold for the kernels . The relations in the first line follows from charge conjugation, whereas those in the second line are obtained by reversing the quark line associated with (but not the one associated with ). The channels on the r.h.s. of (5) are not discussed further in the following. We also recall that kernels for specific flavour transitions are obtained from (4.1) and its analogue for , and that kernels for channels with an initial are equal to the kernels for the charge conjugate channels with an initial .
Useful building blocks for presenting the explicit kernels are the functions
[TABLE]
which are proportional to the LO DGLAP evolution kernels away from the endpoint . The LO kernels and read
[TABLE]
We find it useful to express the NLO kernels in terms of the overcomplete set of momentum fractions and . This makes symmetry properties more transparent and often allows for more compact expressions. Furthermore we will use the notation . Let us first discuss the evolution kernels .
evolution kernels.
Starting with the pure gluon channel, we introduce auxiliary functions
[TABLE]
in terms of which the evolution kernel can be written as
[TABLE]
Here we dropped the subscripts on and for brevity. If we replace and express as a function of only and , we have
[TABLE]
We see in (5) that, apart from the distribution term involving , the kernel is fully symmetric in the momentum fractions of the three final-state gluons. The symmetry under the interchange is trivial, but the symmetry between an observed gluon ( or ) and the unobserved one () is not. We discuss this “active-spectator symmetry” (where we refer to as active partons and to as the spectator) in section 5.2.
For the channel, we define
[TABLE]
and have an evolution kernel
[TABLE]
The evolution kernel for the channel can be expressed by the same functions as
[TABLE]
As we did in the pure gluon channel, we observe active-spectator symmetry for the splitting process . We obtain the kernel from if we drop the distributions terms, omitting and replacing with , and then interchange .
Turning to the channel , we define
[TABLE]
and have an evolution kernel
[TABLE]
In terms of the same functions, the evolution kernel for reads
[TABLE]
This is symmetric in the arguments and , as it must be. Again we observe active-spectator symmetry, this time under the exchange .
For quark-antiquark transitions, we have
[TABLE]
with
[TABLE]
and for the quark valence kernels we have
[TABLE]
with
[TABLE]
The evolution kernels for the remaining channels again obey active-spectator symmetry and are given by
[TABLE]
Terms originating from LO kernels
As specified in (2.3.1), the kernels and can be constructed from and from additional terms that originate from LO kernels. For LO channels, we find that these additional terms have the form
[TABLE]
where or . Note that repeated parton indices on the r.h.s. are not summed over. The flavour diagonal DGLAP kernels in (132) contain distribution terms, which we wish to make explicit. To this end, we write the kernels in the form
[TABLE]
in analogy to (110), where the regular terms read
[TABLE]
We then have convolutions
[TABLE]
where we have omitted the parton labels on for brevity. To obtain these relations, we used equation (B13) in [88] to bring all plus distribution terms into the form .
The convolution terms that appear in the different NLO channels are
[TABLE]
and
[TABLE]
In all these cases, only the flavour non-diagonal DGLAP kernels
[TABLE]
appear, after we used charge conjugation symmetry to replace and .
We note that for all LO and NLO channels , there is exactly one parton combination in each type of convolution term. This will no longer hold for kernels at order , where for given , , and there is more than one possibility for the spectator partons.
In contrast to what we observed for the kernels , active-spectator symmetry does not hold for the LO induced terms just discussed, and therefore this symmetry does not hold for and either. This is most easily seen for . Apart from distribution terms, all functions in the convolutions (2.3) are equal to in that case. The exchange then leaves the first term in (2.3) invariant, whilst interchanging the second and third terms. Because the latter enter with opposite signs in (2.3.1), active-spectator symmetry is broken.
Nonlogarithmic terms.
The full expressions of the kernels , which are not accompanied by logarithms in , are rather lengthy, and we only discuss their general features here. Their full expressions are given in the ancillary files associated with this paper on https://arxiv.org. We note that the functions are not needed for the position space matching kernels , but they do appear in the matching kernels between momentum and position space DPDs.
The general form of the coefficients can be written as
[TABLE]
and
[TABLE]
where the functions with are regular as specified below (110). They are independent of colour factors and of the number of active quark flavours. We see that the colour structure of is the same as for the corresponding kernel .
The regular parts are by far the most lengthy ones, containing rational functions as well as the product of rational functions with logarithms, with products of two logarithms, and with dilogarithms. The logarithms have arguments , , , , , or . The arguments of the dilogarithms can be reduced to the set , and , using the relations given in section II.B of [96]. Note that with these arguments, all logarithms and dilogarithms are real valued. In the denominators of rational functions, the highest powers of and (with ) are the same as in the evolution kernels given earlier.
The plus distribution parts of the coefficients are either zero or simple polynomials in and . They determine the leading behaviour in the limit and will be explicitly given in (5.3.1) below.
The delta distribution parts are in general less lengthy than the regular parts, but they involve rational functions and the product of rational functions with trilogarithms, with products of logarithms and dilogarithms, with dilogarithms, and with products of up to three logarithms. The arguments of all logarithms and polylogarithms are only or . The rational functions often coincide with the corresponding functions in (5), or at least they have the same denominators as these functions.
5.1 Parton number and momentum sum rules
Sum rules for momentum space DPDs at were proposed in [53]. They are among the few theoretical constraints on DPDs that we have. A detailed proof of these sum rules was given in [54], where it was also shown that the renormalisation scale independence of the DPD sum rules implies corresponding sum rules for the evolution kernels. These read
[TABLE]
Note that the distribution terms in the DGLAP evolution kernels on the r.h.s. cancel. In the number sum rule (142), such terms only appear if , in which case the sum of Kronecker deltas gives zero. In the momentum sum rule (143), distribution terms are removed by the prefactor .
For the LO kernels, the sum rules (142) and (143) are readily verified using the relation (37) and the list of possible transitions . At NLO, however, the sum rules provide nontrivial relations between the evolution kernels and the DGLAP splitting functions, which are well known at that order [44, 45, 46, 47, 48, 49, 50, 51, 52]. These relations serve as a valuable cross check of our results.
Number sum rule.
Let us first consider the number sum rule (142) for the different parton combinations. Starting with the ones that involve evolution kernels in LO channels, we have
[TABLE]
where the arguments of the kernels and integration boundaries are as in (142). Here and in the following we omit NLO kernels that vanish identically, such as and . Turning to the NLO channels, we use the notation
[TABLE]
for DGLAP kernels [97] and note that at NLO the relation
[TABLE]
holds for the flavour singlet parts. Taking linear combinations of the number sum rules for the different transitions and and using the symmetry relations (5), we obtain
[TABLE]
We checked explicitly that our results for the kernels fulfil all of the above sum rules.
Momentum sum rule.
For NLO kernels, one finds the momentum sum rules
[TABLE]
all of which are fulfilled for the evolution kernels we have calculated.
5.2 Active-spectator symmetry
In this subsection, we consider kernels away from the kinematic point , so that only real emission graphs contribute and distribution terms do not contribute. We find active-spectator symmetry to hold for the evolution kernels in all parton channels, but not for the any of the kernels that make up the matching coefficient .
The lack of this symmetry for is not surprising, because the renormalisation (2.3) of clearly treats observed and spectator partons in an asymmetric way. Furthermore, already in bare graphs the active partons are singled out because their transverse momenta differ by in the amplitude and its conjugate. This is in fact the only difference between active partons and spectators. To see this, we consider the bare DPDs in light-cone perturbation theory, which can be obtained from covariant perturbation theory by integrating over all minus momenta (see e.g. chapter 7.2.3 of [70] and references therein). Both active and spectator partons are then on their mass shell by construction. As shown in section 5 of [54], for the numerator factors in light-cone perturbation theory are also identical for active and spectator partons in bare DPDs, provided of course that one considers unpolarised partons.
To understand why active-spectator symmetry holds for , we note that for these evolution kernels are associated with the ultraviolet divergences of two-loop graphs for the splitting process , where is the spectator parton. More precisely, these kernels arise from kinematic configurations in which all three final state partons have large transverse momenta. The finite value of can be neglected in that region, so that one should have active-spectator symmetry according to our arguments in the previous paragraph. Note that divergences from configurations in which only two of the three final state partons have transverse momenta in the ultraviolet are not associated with , but with the lower-order splitting kernels or .
Let us give another reason why the finite transverse momentum in our calculation should play no role for the kernels . We obtained from the computation of the splitting kernels , where plays the role of a hard scale and cannot be set to zero. On the other hand, the evolution kernels appear in the evolution equation of the DPDs , which is valid at . One could hence compute from graphs with , provided that one has a dimensionful infrared regulator for separating infrared and ultraviolet poles in dimensional regularisation.
5.3 Kinematic limits
In this subsection, we analyse the behaviour of the kernels and in various kinematic limits.
5.3.1 Threshold limit: large
To begin with, we consider the convolution (5) at the parton-level threshold , i.e. . Writing out the plus-distribution, we readily obtain
[TABLE]
A logarithm in is thus generated by the plus distribution term of the kernel , which hence dominates the threshold behaviour of the convolution. Since the kernels do not contain distributions with , no higher powers of are generated in the convolution (5). This situation is analogous to the convolution of ordinary DGLAP splitting functions with PDFs, where no such distributions appear even at NNLO [98, 99].
The plus distribution parts of the LO channel kernels are given by
[TABLE]
and
[TABLE]
where we omitted coefficients and that are exactly zero.
The kernels and for NLO channels have no plus distribution terms. Their leading threshold behaviour is with , which gives a convolution of order . Specifically, we have
[TABLE]
and
[TABLE]
5.3.2 Small
We now turn to the case in which both momentum fractions and are small; this is a typical situation at high collision energy when two systems of moderately large invariant mass are produced. In the limit the convolution integral (5) includes the region in the kernel , where is defined in terms of the external momentum fractions and by (113).
Let us first recall the analogous situation for the one-dimensional Mellin convolution of an LO DGLAP kernel and a PDF,
[TABLE]
If and for , then the region in (154) gives a behaviour proportional to
[TABLE]
where we obtained the leading logarithmic behaviour in on the r.h.s. by extending the integration to the full range . This is nothing but the generation of small- logarithms in a PDF: starting with powers of , the convolution with a kernel proportional to results in an additional power of that logarithm.
In full analogy one obtains a small- logarithm in the two-variable convolution (5) if has the same small- behaviour as above and
[TABLE]
The region then contributes to as
[TABLE]
where on the r.h.s. we extracted the leading logarithmic behaviour by extending the integration to the full range . If grows less fast than at small , then the region is no longer dominant in the convolution integral, and no small- logarithm is built up by the integration over . Notice also that the dependence of the kernel in (156) directly determines the dependence of the convolution (157).
For the evolution kernels , we find the following leading behaviour at small :
[TABLE]
For the kernels without a singularity, we find
[TABLE]
For the coefficients going like at small , we find
[TABLE]
and
[TABLE]
while for the coefficients with a less singular behaviour we have
[TABLE]
For both types of kernels, the channels with a behaviour are , , , and . We checked that the graphs giving rise to such a behaviour are those in which — on both sides of the final-state cut — the spectator parton is emitted from a three-particle vertex at which a slow gluon is emitted, where “slow” is relative to the incoming parton at that vertex. Examples for such graphs are figures 1(a) to 1(d) and figures 2(g) to 2(j).
Notice finally the Casimir scaling between and , and between and . The small limits of , , and (but not of ) in these pairs of channels are simply related by a proportionality factor .
Triple Regge limit.
The results just presented also allow us to analyse the “triple Regge limit” . To this end, we simply take the limit in (157). Restricting our attention to the kernels with a behaviour, we have
[TABLE]
and
[TABLE]
and
[TABLE]
We see that the leading singular behaviour of all these kernels involves at most two powers of . An exception are the kernels for the emission of two gluons. In this case we have a power-law behaviour like of the kernels, which results in a power-law behaviour like of the convolution (157). One finds that this behaviour comes from the terms and and corresponds to graphs with UD topology as shown in figures 1(b) and 2(h). In these graphs there are two consecutive three-point vertices at which an observed slow gluon is radiated from a parton carrying the full or almost the full initial plus-momentum. The appearance of a power behaviour only in but not in or means that this behaviour goes along with a double logarithm . This corresponds to strong ordering in the transverse momenta at the two consecutive splitting vertices, which gives two logarithmic transverse-momentum integrals.
Analogous expressions for limit are easy to obtain, because the channels with a behaviour have kernels that are symmetric in and .
5.3.3 Small or
We now consider the limit in which whilst is not. This is relevant for DPS processes in which one of the two hard systems is very heavy, whereas the other one is light compared with the total collision energy. The integration in the convolution (5) can in this case not reach the region . Hence, the integration cannot build up a small- logarithm, and the limit of the convolution is simply obtained from the limit in the kernel . Note that, unlike in the case , the distribution parts and of the kernel can also contribute to the small behaviour.
In the following, we give the behaviour of all kernels for the case in which whilst is unconstrained. This covers the limit just discussed, and it will also allow us to consider the nested limit later on. An analogous discussion holds of course for the limit at generic values of . For small and generic , one finds for the evolution kernels
[TABLE]
while the ones without a singularity go like
[TABLE]
The coefficients are in this limit given by
[TABLE]
and
[TABLE]
The limiting behaviour of the subleading kernels is given by
[TABLE]
We see that the channels with a singularity are those in which the parton with momentum fraction is a gluon. There are hence graphs in which that slow gluon is radiated from a parton with momentum fraction much larger than , i.e. one has a vertex with the emission of a “slow” gluon in the same sense as specified after (5.3.2).
The preceding expressions are obtained by approximating the kernels for . One may subsequently take the limit , which becomes relevant in the convolution if . Comparing the results of this procedure with those in (5.3.2) to (5.3.2), we find that in general the limits and do not commute. The only case when they do commute is if a kernel has the maximally singular behaviour in both limits. This holds for the kernels and we already discussed after (5.3.2). In all other cases, the behaviour of the convolution in the triple Regge limit depends on the direction in which this limit is approached.
5.4 Numerical illustration
A systematic investigation of the numerical impact of the NLO kernels we have computed is a substantial task, far beyond the scope of this paper. Let us, however, illustrate the impact of the NLO corrections to splitting in a simple setting. We compute the perturbative splitting contribution to the two-gluon DPD at NLO,666Note that the label “NLO” on refers to the sum of and terms.
[TABLE]
with
[TABLE]
where and the parton densities on the r.h.s. are evolved at NLO. We set , so that only the part of the second-order splitting kernels contributes. We take for the transverse distance between the partons, which corresponds to . We use the CT14 PDF set [100] at NLO, with an associated value of for the strong coupling at the scale considered here.
In figure 7 we show along with its first-order contribution . We also show the leading-order expression of the two-gluon DPD, which has the same form as in (172) but with and the gluon distribution evaluated at LO. The latter are again taken from the CT14 analysis, with an LO coupling of .777We use the sets CT14llo and CT14nlo from LHAPDF [101] for and PDFs at LO and NLO, respectively. In figure 8 we show the ratios
[TABLE]
The first ratio quantifies the relative change when going from the LO to the NLO approximation of the DPD, whereas the second ratio indicates the relative importance of the kernels in the NLO results. Notice that in figure 8 goes up to the kinematic boundary .
We see that the difference between at LO and NLO can be appreciable, especially for small momentum fractions and when is close to its kinematic boundary . Typically, there is a considerable difference between and , which simply reflects the difference between at LO and NLO. That this difference is large for the gluon density at moderately large scales is in fact well known and not specific to the CT14 parton distributions. Compared to this effect, the impact of the term is smaller, but nevertheless its size amounts to to of in a wide kinematic regime. We also produced analogous plots with the LO and NLO sets of MMHT 2014 [102] and of NNPDF3.0 [103]. For up to 0.1, they look quite similar to those for CT14. For larger momentum fractions the differences between the sets become more pronounced, due to larger differences between the respective gluon distributions.
In figure 8 we also show the ratio
[TABLE]
with
[TABLE]
which quantifies the contribution of the splitting processes and to . Comparing with , we see that the contribution from is small compared with the one from , except when approaches its upper boundary at fixed . Interestingly, both for and , there is a value of just below at which the contributions from and cancel each other, whilst each of them separately is of order .
5.5 Comparison to other results in the literature
Various papers in the literature contain results for QCD splitting functions at that are (at least) differential in the momentum fractions of two of the produced partons. However, these functions have been derived in somewhat different contexts for other observables. One might wonder whether these objects and the evolution kernels we have computed here are the same. In general, one should not be surprised if there is a difference at even if the functions agree at — recall that in the single inclusive case, the DGLAP splitting functions appropriate to PDFs (spacelike case) and fragmentation functions (timelike case) differ at [48, 49]. It is possible to connect the two beyond LO [104, 105, 106, 107, 48, 49, 50, 108, 109, 110, 111], but such connections are nontrivial and go beyond simple equality.
In [112, 113], the corrections to splitting functions relevant to the evolution of so-called fracture functions were computed for most partonic channels (excluding, in particular, ). These fracture functions appear in semi-inclusive DIS when we measure the energy of a hadron in the forward region. They have two momentum fraction arguments, one for the parton entering the hard scattering process, and the other for the measured forward hadron. The associated splitting functions thus involve the measurement of the momentum fraction of one spacelike parton associated with a PDF operator, and that of a timelike parton that ultimately yields the final-state hadron. By contrast, for our DPD splitting functions, we measure the momentum fractions of two spacelike partons.
In [114] it was proposed that the NLO splitting functions in the evolution of fracture functions are actually equal to the evolution kernels for DPDs, after a simple transformation of the momentum fraction arguments to take into account the fact that one of the momentum fraction arguments in the splitting functions of [112, 113] is not the light-cone momentum of the final parton over that of the initial parton, but that of the final parton over that of the other final parton. We find that this equality does not hold, and that in general the transformed splitting functions of [112, 113] do not have the symmetry that our kernels have (and must have) under the interchange of partons and . For example, in our results is symmetric under , which is due to the equivalent status of partons and together with charge conjugation invariance. By contrast, the same is not true for the corresponding function in [112, 113].
References [115] and [116] give splitting functions at for the decay of a timelike parton, where one measures the momentum fractions of two of the partons produced. These functions are referred to as “two-body inclusive decay probabilities” in [115, 116]. They appear in the evolution of two-particle fragmentation functions, which in turn are relevant when one measures the energies of two hadrons in an individual jet [117]. We find that in general these functions also differ from the DPD kernels . Exceptions to this are the valence-type kernels and , both of which agree between the spacelike and timelike cases. This is likely due to the fact that for each of these kernels there is only one graph in light-cone gauge, see figure 2(m) and 2(o), which does not have any infrared or ultraviolet sub-divergences. This finding is reminiscent of the equality of the single inclusive splitting function between the spacelike and timelike cases [48].
Finally, we comment on the tree-level “triple collinear” splitting functions computed in [118, 119]. These quantities are completely unintegrated, depending on the total invariant mass , the invariant masses of parton pairs , as well as the momentum fractions of all partons. They are thus analogous to the initial real emission graphs we compute using FORM or FeynCalc. However, in our calculation we require a momentum shift between amplitude and conjugate amplitude for the two observed partons, whereas the functions computed in [118, 119] have no such shift. Thus we could not have used these functions as the starting point of our calculation for the real emission diagrams. Similar statements apply regarding our virtual diagrams and the one-loop double collinear splitting functions [120, 121, 122, 123, 124, 125, 126], which describe the exclusive splitting process at the one-virtual-loop level. It should be noted that [120, 121, 122, 123, 124, 125, 126] also give results for one-loop double collinear “splitting amplitudes”, which describe the one-loop splitting process on one side of the cut. We could have potentially used these amplitudes in our calculation. It is worth mentioning that the tree-level triple collinear splitting functions and one-loop double collinear splitting functions can be used as the starting point in many NNLO beam and jet function computations, which do not require a momentum shift between amplitude and conjugate, as was pointed out in [127].
5.6 The SUSY Relation
There is an interesting relation between the gluon- and quark-initiated single-inclusive splitting functions at if we set and . Namely, if we define
[TABLE]
then we have [128]. At NLO we have [49], but the expression for is considerably shorter than those for the individual splitting functions.
This relation holds at LO because when one modifies the colour factors as specified, the theory becomes equivalent to an supersymmetric theory, with the quark playing the role of the gluino [129]. We have beyond LO because the conventional methods of dimensional regularisation and renormalisation (which are used in [49]) do not preserve supersymmetry. Since the relation is only violated by the regulator, the expressions for at higher orders are still comparatively small. If one instead uses a regularisation method that preserves supersymmetry, such as dimensional reduction [130], then holds beyond LO, as was shown for the NLO case in [131, 132, 133, 134, 108].
It is interesting to see if a corresponding relation is “approximately” observed by our evolution kernels, where we recall that these have also been obtained using dimensional regularisation and renormalisation. The straightforward generalisation of the difference (177) to this case is
[TABLE]
The relation between and readily implies that . With the NLO kernels given earlier in this section, we find
[TABLE]
where it is understood that the plus prescription in the first term is dropped when is replaced by or . This is indeed a much simpler expression than those for the individual evolution kernels. Note in particular that for , the logarithmic terms completely cancel in , and only the single logarithmic terms survive at . Presumably we would obtain if we were to convert our results to the dimensional reduction scheme, but to verify this goes beyond the scope of the present study.
6 Conclusion
We computed at two-loop order the perturbative matching kernels between position space DPDs and PDFs , as well as the matching kernels between momentum space DPDs and PDFs . The computation of is more involved than the one of because it requires results for bare graphs at one higher order in . We showed that at two-loop accuracy one can readily extract from the kernels , the evolution kernels in the inhomogeneous term of the evolution equation for , as well as the matching coefficients for computing from . We obtained results for all possible partonic channels, i.e. all possible combinations of , and , while limiting ourselves to unpolarised, colour singlet DPDs. These quantities are needed to compute DPS cross sections at next-to-leading order.
We performed the calculation using both Feynman and light-cone gauge, finding agreement between the results in the two gauges. Loop integrals were performed using the method of differential equations, together with boundary conditions computed using the method of regions and standard integration using Feynman parameters. Each master integral was checked numerically at values of .
We extracted the behaviour of the matching coefficients in the threshold limit and in different limits with small momentum fractions, namely for , for , and for . We verified that our two-loop results for obey the number and momentum sum rules given in [54]; this can be viewed as a cross check of our results, or as an explicit verification of these sum rules at the two-loop level. In [114] a relation was proposed between the kernels and the splitting functions appearing in the evolution of fracture functions. We found that this relation is not fulfilled at two-loop order.
It is quite straightforward to extend the results of the present paper to polarised DPDs. Similar methods as shown here can be used to compute the matching coefficients for colour interference. In this case, one needs to explicitly deal with rapidity divergences and soft factors, as briefly described at the end of section 4.2.1. Work in this direction is underway. On the phenomenological side, the numerical example in section 5.4 shows that the difference between two-gluon DPDs computed at leading and at next-to-leading order can be appreciable, in particular at small momentum fractions and for close to . It will be interesting to investigate in a more comprehensive way how the two-loop corrections computed in this paper affect DPDs, double parton luminosities, and ultimately the cross sections of processes that are sensitive to DPS.
Acknowledgements
We gratefully acknowledge discussions with Johannes Michel. Two of us (MD and JRG) thank the Erwin Schrödinger International Institute for Mathematics and Physics (ESI) for hospitality during the Parton Showers, Event Generators and Resummation Workshop (PSR19), when portions of this work were completed. The figures in this work were produced with JaxoDraw [135, 136].
We would like to dedicate this work to the memory of James Stirling.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. Kulesza and W. J. Stirling, Like sign W 𝑊 W boson production at the LHC as a probe of double parton scattering , Phys. Lett. B 475 , 168 (2000), 10.1016/S 0370-2693(99)01512-9 , hep-ph/9912232 . · doi ↗
- 2[2] J. R. Gaunt, C.-H. Kom, A. Kulesza and W. J. Stirling, Same-sign W 𝑊 W pair production as a probe of double parton scattering at the LHC , Eur. Phys. J. C 69 , 53 (2010), 10.1140/epjc/s 10052-010-1362-y , 1003.3953 . · doi ↗
- 3[3] CMS Collaboration, Double Parton Scattering cross section limit from same-sign W bosons pair production in di-muon final state at LHC , CMS-PAS-FSQ-13-001 (2015).
- 4[4] CMS Collaboration, Measurement of double parton scattering in same-sign WW production in p-p collisions at s = 13 Te V 𝑠 13 Te V \sqrt{s}=13\leavevmode\nobreak\ \mathrm{Te V} with the CMS experiment , CMS-PAS-FSQ-16-009 (2017).
- 5[5] F. A. Ceccopieri, M. Rinaldi and S. Scopetta, Parton correlations in same-sign W 𝑊 W pair production via double parton scattering at the LHC , Phys. Rev. D 95 (11), 114030 (2017), 10.1103/Phys Rev D.95.114030 , 1702.05363 . · doi ↗
- 6[6] S. Cotogno, T. Kasemets and M. Myska, A spin on same-sign W-boson pair production (2018), 1809.09024 .
- 7[7] R. Aaij et al. , Measurement of the J/ ψ 𝜓 \psi pair production cross-section in pp collisions at s = 13 𝑠 13 \sqrt{s}=13 Te V , JHEP 06 , 047 (2017), 10.1007/JHEP 06(2017)047 , [Erratum: JHEP 10,068(2017)], 1612.07451 . · doi ↗
- 8[8] V. Khachatryan et al. , Search for new physics in same-sign dilepton events in proton–proton collisions at s = 13 Te V 𝑠 13 Te V \sqrt{s}=13\,\text{Te V} , Eur. Phys. J. C 76 (8), 439 (2016), 10.1140/epjc/s 10052-016-4261-z , 1605.03171 . · doi ↗