Two real parton contributions to non-singlet kernels for exclusive QCD DGLAP evolution

Results for the two real parton differential distributions needed for implementing a next-to-leading order (NLO) parton shower Monte Carlo are presented. They are also integrated over the phase space in order to provide solid numerical control of the MC codes and for the discussion of the differences between the standard $\bar{MS}$ factorization and Monte Carlo implementation at the level of inclusive NLO evolution kernels. Presented results cover the class of non-singlet diagrams entering into NLO kernels. The classic work of Curci-Furmanski-Pertonzio was used as a guide in the calculations.


Introduction
With the start of the operation of the Large Hadron Collider (LHC) in CERN and the progress in the analysis of growing data samples for many scattering processes, mastering the precise evaluation of strong interactions effects within perturbative Quantum Chromodynamics (QCD) [1][2][3] will become quickly more and more important. QCD is a well established theory -testing the validity of QCD is not an open issue any more. Its principal role in the LHC data analysis will be providing precise predictions for rates and distributions of quarks, gluons and hopefully other newly found particles carrying colour charge, being part of either signal process or background. The most important theoretical tools in perturbative QCD (pQCD) calculations, apart from Feynman diagrams, renormalization, etc. are the so called factorization theorems, see for instance [4][5][6], which allow to describe any scattering process with a single large mass or transverse momentum scale µ F (enforcing short distance interaction), in terms of the on-shell hard process matrix element (ME) squared and convoluted with the ladder parts. The hard process is calculated up to a fixed perturbative order. The ladder parts are defined for each colored energetic parton entering (exiting) the hard process. The initial state ladders are conveniently encapsulated in the inclusive parton distribution functions, PDFs.
The logarithmic dependence of the parton distribution functions (PDFs) on the large scale µ F is described as a DGLAP [7] evolution of the PDFs. This evolution was studied for the inclusive PDFs up to NLO level in the early 80's, see refs. [8,9], and was recently established at the NNLO level [10,11].
Instead of being encapsulated in the inclusive PDF, the multi-parton emission process of the initial state ladder, can be modelled using direct stochastic simulation in terms of four-momenta and other quantum numbers, within the Monte Carlo (MC) parton shower. Here, the baseline works have been done in mid-80's, see refs. [12,13], where the leading order (LO) ladder was implemented in parton shower MC (PSMC) programs. Standard LO level PSMCs implement also the process of hadronization of the light quarks and gluons into hadrons and play an important role in the software for all collider experiments because of that.
Fulfilling the challenging requirements on the quality and precision of the pQCD calculations, needed for the experimental data analysis at the LHC, enforces for the first time an urgent solution of the problem of upgrading exclusive PSMC to the complete NLO level, the same level, which was reached for inclusive PDFs two decades ago. This task is highly nontrivial mainly because the classic factorization theorems [4][5][6] were never designed for the exclusive MC implementation, but rather for defining inclusive PDFs and performing fixed order calculations for the hard process, convoluted with these PDFs.
Let us comment briefly on the longer term physics impact of the present work. Remembering that QCD is not any more a new theory, the main impact of this work will be the improvement of calculations of QCD effects, for hadron collider experiments like LHC, with the aim of improving chances of discovering directly or indirectly New Physics and/or better measurement of the Electroweak Standard Model parameters, especially when high statistics, high precision data are accumulated.
More precisely, this work elaborates on pQCD effects in the initial state, which from the perspective of the LHC experiments influence mainly: (A) overall normalization of the hard processes through parton luminosities, (B) distributions of transverse momenta (k T ) of incoming partons deforming many other important distributions in any hard process, including searches for supersymmetry, etc., (C) and provide one or more jets accompanying hard process.
The longstanding problem in the data analysis is that the above three classes of important QCD effects are addressed by three separate theoretical pQCD calculational tools, based on different incompatible perturbative techniques: (A) strictly collinear NLO PDFs [7], (B) semi-inclusive schemes of infinite order soft gluon k T resummation [5,14,15] or, alternatively, LO parton shower Monte Carlos [12,16] (C) finite order NLO (NNLO) calculations [17,18] sometimes combined with a LO parton shower MC [19] or collinear PDFs.
In the analysis of experimental data combining these (and other) techniques is not only inconvenient, but also is a serious source of irreducible theoretical (and experimental) uncertainties. This old and well known problem becomes more severe with the increasing precision of the collider data, and will inevitably plague high statistics, high precision LHC data. The ultimate aim of the present work is to provide a basis for designing a single Monte Carlo program addressing all three classes of the above pQCD effects at once, within the same consistent pQCD theoretical framework. However, for this to be realized one has to start with solving one basic difficult problem -the upgrade of the initial state parton shower MC to at least the same level as standard PDFs, that is to the NLO level, in the fully exclusive manner. The present work provides essential building blocks (exclusive kernels) for this critical extension of the parton shower, and analyzes factorization scheme differences with respect to the standard CFP M S scheme. While this work focuses on the implementation of NLO corrections to the initial state ladder parts (parton showers), the hard process part at NLO within the same scheme is discussed in ref. [20] 1 .
The critical problems to be solved on the way to NLO PSMC are the following: 1. Violation of four momentum conservation. In the standard collinear factorization four momentum conservation is broken both between the hard process and the ladder as well as between the ladder segments (2PI kernels). The source of this non-conservation in the standard collinear factorization is the introduction of the projection operators, P kin in ref. [4], operator P in ref. [9], or operator Z in ref. [24].
These operators are absent in the first step of the separation of the collinear singularities into the ladder parts -they are introduced later on in order to (i) conveniently isolate the lightcone variable integration out of the phase space for analytical integration and (ii) facilitate the order-by-order pQCD calculations separately for the hard process ME and the ladder elements (kernels). This non-conservation happens in the transverse momenta, which are anyway treated inclusively (integrated over), hence this bad feature of collinear factorization goes almost unnoticed. 2 In the traditional LO PSMCs the above non-conservation is repaired "by hand" [12,13], but in such a way that the NLO effects induced by this reparation are analytically almost uncontrollable. This is not a problem, unless one attempts to complete NLO in the ladder, or to combine an NLO hard process ME with a LO PSMC [19]. A systematic solution of the above problem must involve replacing the projection operator P kin by a more sophisticated operation involving a special parametrization of the entire phase space for the hard process and the ladders. An explicit example of such a solution for the W/Z production process (Drell-Yan process) can be found in ref. [20].
2. Separation of singularities before phase-space integrations. In the collinear factorization, separation of the LO singular contributions in the form of the leading logarithm ln Q 2 m 2 or 1 ε poles of dimensional regularization can be done only after the phase space integration. In order to construct an efficient Monte Carlo algorithm this separation has to be done at the very beginning, at the integrand level, before the phase space integration. This requires going beyond the inherently inclusive approach of collinear factorization. The effort of getting the NLO prediction for the semi-inclusive distribution in the phase space of the hard process (like k T and rapidity distributions in W/Z production) started already quite early, see for instance [17,25]. More recently, Monte Carlo tools combining NLO ME for the hard process with the LO PSMC were developed in [19] and [26]. Studies on redefining PDFs in a partly exclusive form beyond LO (unintegrated PDFs) [5], or in exclusive form (fully unintegrated PDFs) [27] are also pursued.
3. Negative "probability distributions". The NLO corrections in collinear factorization are negative in some regions of phase space and therefore cannot be generated directly in the Monte Carlo, if we insist on the realistic simulation using positiveweight MC events. 3 The reasons for non-positiveness of the NLO corrections is well understood. For instance, in the physical gauge positive squares of the Feynman diagrams are typically more divergent and form the LO approximation, while nonpositive interferences are collected in NLO corrections. The use of kinematic projection operators in factorization and related subtractions are another sources. The factorization procedure collects all these non-positive corrections into separate objects and the integration over the phase space is done for each of them separately. Consequently, part of the factorization procedure has to be reversed in order to recombine non-positive NLO distributions with the positive LO distributions, before the MC algorithm is designed. The above defactorization procedure has to be outlined. An example proposal relying on Bose-Einstein symmetrization was formulated in ref. [28] and an even more promissing one is described in ref. [29] (similar to that in ref. [30]).

4.
Lack of the published exclusive NLO distributions. We have not found exclusive NLO distributions forming the NLO corrections in the ladder part in literature -all published results in the context of the NLO calculations of kernels for evolution of PDFs are integrated over the phase space. The main objective of this paper is to provide such distributions for the non-singlet case.
5. Inclusive treatment of multigluon soft limit. One of the critical issues in any type of collinear factorization is the behavior of many-gluon distributions in the soft limit. In the inclusive approach it is enough to know that the cancellations between real and virtual soft contributions allow us in principle to neglect entire classes of diagrams and/or divergent contributions -they sum up to zero. In the exclusive MC approach these contributions/singularities, instead of being dropped out, have to be modelled precisely to infinite order. The above soft gluon behavior is rather well known, albeit more complicated in QCD than in QED due to the (non-abelian) triple gluon vertex -the eikonal limit and angular ordering are known to govern it [31]. In ref. [32] a detailed analysis of the soft limit for the C 2 F and C F C A contributions for non-singlet evolution kernels was performed and the well known angular ordering was exposed, both numerically and analytically, for the first time using exact double gluon distributions presented explicitly. 6. Inappropriateness of M S factorization scheme. The issue of factorization scheme dependence has to be revisited and one has to decide about the best factorization scheme (FS) for the PSMC implementation. It should not be taken for granted that it will be the M S scheme, which is presently the "industry standard" for the inclusive PDFs and in pQCD calculations for the hard process. 4 In the view of the above the following question emerges: Can one stay strictly within the M S scheme for the MC implementation of the NLO ladders combined with the NLO hard process ME? In the present work we will addres this question partly, for the MC modelling of the ladder part. It is studied for the hard process ME in ref. [20].
Let us indicate that our answer will be negative: FS of the NLO MC has to differ from M S scheme for several reasons. Some of them will be discussed in detail in this paper. The key element defining FS are the so-called soft collinear counterterms. 5 In MC we choose them to be identical with the distribution used in the LO MC, e.g. the LO level PSMC is constructed by means of iterating soft collinear counterterms. This choice on one hand will simplify NLO corrections in the MC, but it will also cause departure from the M S FS. The second source of discrepancy will be the fact that in the MC, the factorization scale will be identical to a well defined kinematic variable Q F of the LO MC, replacing the formal parameter µ F of the M S FS. The logarithm of Q F is then used in the MC as an ordered evolution time variable. Q F will typically be maximum transverse momentum, maximum angle, or virtuality of the emitted particles. For the gluonstrahlung, the choice of the maximum transverse momentum results in the MC FS very close to M S FS. 6 However, the soft gluon eikonal limit for contributions like gluon pair production or quark-gluon transitions, dictates the use of the variable related to the maximum angle (rapidity) as a factorization scale variable in the MC. Third reason is the presence in the M S FS of certain artifacts of dimensional regularization which cannot be implemented in the MC in four dimensions.
7. Constrained evolution. Constrains on the momenta and other quantum numbers imposed by the hard process on the initial state parton shower (ladder), especially important in the presence of the narrow resonances, must be taken into account by the PSMC. The clever and powerfull MC technique of backward evolution of ref. [12] is good for the LO PSMC, but for the NLO case one may need something more sophisticated. The dependence of the ladder part (PDF) on the factorization scale is described in pQCD by the integro-differential DGLAP [7] evolution equation. Its 4 Since then the FS dependence is reduced to the discussion of the residual dependence on the factorization scale µF of M S scheme due to higher orders. 5 We refer to the exclusive version of soft collinear counterterm, which is the distribution within the 1-particle phase space encapsulating collinear (and soft) singularity -not its integral as in inclusive FS. 6 This fact is well known, see also analysis of ref. [33].
solution has the form of a time ordered (T.O.) exponential with the ordering in the logarithm of the factorization scale. 7 This T.O. exponential can also be obtained directly from the ladder Feynman diagrams. In the Monte Carlo T.O. exponential is conveniently modelled using a Markovian MC algorithm. 8 However, a Markovian MC algorithm would be highly inefficient for modelling the initial state radiation (ISR) ladder, when hard process selects very narrow range of energies and flavors of the partons incoming into hard process (with narrow W/Z resonances). Here, the backward evolution MC algorithm is a standard solution [12] -used in all standard LO PSMCs. It is to be seen whether backward evolution MC is upgradable to NLO PSMC. In the meantime, the constrained MC algorithm of refs. [35,36] offers an interesting alternative. While the backward evolution MC needs pretabulated solutions (PDFs) of the evolution equations, which have to be prepared beforehand using non-MC auxiliary codes, the constrained MC performs evolution on its own.
In this article we mainly addres point 4 in the above list of problems. However also point 6 is discussed in some detail.
Let us now establish precise terminology concerning diagrams and phase space integration. We will elaborate on the diagrams contributing to non-singlet evolution kernels at the NLO level.
Generally, in the calculation of the exclusive/inclusive evolution kernels in this work we will take paper of Curci Furmanski and Petronzio [9] (CFP) as the starting point and as the reference in all NLO calculations. For the non-singlet part of the QCD DGLAP evolution we will calculate (or recalculate) both the standard inclusive NLO kernels and the new exclusive (unintegrated) ones, which are needed for constructing NLO PSMC. This work prepares building blocks for NLO PSMC, whereas the actual MC algorithm and all issues related to the factorization scheme used in NLO PSMC will be discussed in a separate paper [20].
Our main object of interest in the present work are the diagrams depicted in Fig. 1, with two emitted on-shell quarks and/or gluons, that is diagrams with two cut lines. We will call them 2-real, or shortly 2R, contributions. The other diagrams with 1-real and 1virtual, nicknamed 1V1R, will be also partly considered. Diagrams with 2 virtual (2V) will not be discussed, because they will be treated in the same way as in CFP (deduced from the baryon number conservation rule, which we keep in the Monte Carlo by construction).
The NLO 2PI diagrams feature amplitude-squares, depicted in the upper row in Fig. 1: double gluon emission 1(a), gluon pair production 1(b) and fermion pair production 1(c) as well as interferences, displayed in the lower row in Fig. 1. Interference diagrams enter into the MC kernels for the first time at NLO and their correct incorporation in the MC requires more care. They can potentially make negative contributions to the kernel, spoiling the MC weight. It turns out that the interference diagrams have different statuses in the MC, in the following we explain the reasons for it. Figure 1. Real contributions to NLO non-singlet DGLAP kernel.
Let us make a simple observation that interferences 1(d) and 1(e) together with squares 1(a) and 1(b) form the full amplitude square 2 .
Positiveness of this amplitude square will cause the MC weight introducing these two interferences to be positive. Both interferences 1(e) and 1(d) can be implemented in the MC in the 2R group of diagrams.
The interference diagram 1(f) originates from the following amplitude-squared where qq pair production 1(c) is already included in the non-singlet class, while qq transitions amplitude (squared) belongs to the singlet kernel. Apparently, diagram 1(f) corrects quark-gluon transitions absent in non-singlet kernel and must be included in the MC together with the singlet contributions. Similarly, the diagram 1(g) is the interference part in the amplitude square 2 , where both amplitudes representing qq transitions essentially belong to the singlet NLO kernel.
In the case of both interferences related to quark-gluon transitions the MC weight is not protected by the Schwartz inequality due to the absence of amplitude squares in the non-singlet kernels. The above problem will naturally disappear when singlet diagrams are included in the game, and we need some temporary fix, while staying in the non-singlet class. In order to keep the fermion number conservation we keep these contributions in the kernel, but add them to the 1R1V class and treat inclusively. In the following such a combination of the 1R1V and integrated subclass of 2R interferences we will call collectively unresolved contribution.
Another important point concerns internal singularities of Feynman diagrams. They are present in graphs of Fig. 1(a), 1(b) and 1(c). The double bremsstrahlung diagram Fig. 1(a) does not enter into NLO kernel as a whole 9 but only what remains after subtracting soft collinear counterterm of the CFP factorization. On the other hand, diagram of Fig. 1(b) features internal collinear singularity cancelled by the corresponding virtual (gluon self-energy) diagram. These diagrams may enter into the unresolved part in the MC, or may be modelled in an exclusive manner. In the latter case diagram 1(b) requires the construction of a dedicated soft-collinear counterterm which encapsulates the above internal singularities and is instrumental in the MC construction. Such a counterterm will be defined in Section 3.6. Diagram 1(c) can be also modelled in an exclusive way and also needs a dedicated counterterm, see Section 4.3.
Before that, Section 2 presents an overview of calculations illustrated by the example of the LO kernel. Notation and methodology of extracting kernels will be introduced using this simple example. The integration procedure for NLO kernels is reviewed in Sections 3 and 4. The contributions from each appropriate Feynman diagram in the axial gauge are calculated in integrated and unintegrated form needed for the MC. Analytical integration for control will also be done. Discussion of the collinear and soft singularity structure will have the highest priority.
Sections 3 and 4 summarize results of analytical integrations. In contrast to the results presented in [9], 2R contributions will be shown separately instead of the sum of 2R and 1R1V. Section 5 provides final discussions and conclusions.

Leading order example
Using the LO diagram of Fig. 2 we are going to introduce some basic notation, which will also be useful in the NLO calculations. In addition we will show correspondence between elements of the CFP [9] scheme using dimensional regularization (n = 4 + 2 ) and calculations in the Monte Carlo (n = 4) in this simple example.
On one hand, the non-singlet bare PDF of CFP [9], from which the DGLAP kernel is extracted, reads as follows It is not two-particle-irreducible (2PI) where PP is the pole part operator, C F is the color factor, Z F = 1 + Z [1] + . . . is the fermion wave function renormalization factor. Let us explain step by step the other elements in the above formula. One particle phase space of emitted on-shell gluon (cut line), in n = 4 + 2 dimensions, reads: (2. 2) The emitted gluon 4-momentum k µ is parametrized using Sudakov variables where p is the momentum of the incoming parton and the lightlike vector n defines axial gauge. The conditions . Non-abelian coherence effects in the soft limit beyond LO are easier to handle if we use the variable instead of transverse momentum k. Its modulus, |a| ≡ a, will be used to define the factorization scale in the MC -we will refer to it as an angular scale. 10 Phase space in eq. (2.1) requires to be closed from the above, at least temporarily. In the CFP scheme this closure plays a marginal role, as Γ(x, ) consists of pure poles. The upper limit merely influences intermediate results, through the parametrization of the phase space, before taking PP and executing all kind of internal infrared (IR) cancellations. On the contrary in the Monte Carlo scheme the variable defining the upper limit of the phase space plays an important role of the factorization scale. Its logarithm is the evolution time variable in the MC algorithm/code. The most popular choices are: maximum angular scale (angular ordering), and maximum transverse momentum of all real emitted partons. Virtuality of the emitter line in the ladder or maximum k − = k 0 − k 3 are the other valid but less attractive choices.
Setting the upper phase space limit for the angular variable, |a| ≤ Q, we obtain: 11 (2.6) In the above Z The evolution kernel is defined as twice the residue of Γ at = 0: How does the above compare with the Monte Carlo? In the Monte Carlo the same integral taken in the limits q 0 < a < Q in n = 4 looks simpler: is the Sudakov formfactor. As we see, the same LO kernel is now the coefficient in front of the collinear logarithm: Apparently, 1/ pole of CFP translates into ln Q q 0 : The relation between PDFs and evolution kernels in CFP and Monte Carlo factorization schemes is, however, more complicated beyond LO. This is discussed in more detail in ref. [37], see also refs. [20,33]. Generally, all differences between MC and CFP factorization schemes will be traced back to diagrams with subtractions, or with internal collinear singularity cancellations.

2R contributions to non-singlet NLO kernels
In the following we shall calculate bare PDFs and the resulting inclusive evolution kernels from 2-real phase space of the Feynman diagrams contributing to the non-singlet NLO DGLAP kernels in QCD. Let us stress that our real aim are the distributions over the 2R phase space integration. Analytical integrations will be performed mainly as a crosscheck with the know results and for testing parts of the Monte Carlo code. We shall start with explaining notation and 2R phase space parametrization used in the calculations. Differential distributions and 2R phase space integrals will be listed for each Feynman diagram separately.

Kinematics
Sudakov parametrization is introduced for both emitted partons: where p µ is the momentum of the incoming quark (p 2 = 0) and lightlike n µ is the axial gauge vector. Real on-shell emitted gluon or quark has 4-momentum k µ i , and q µ i denotes 4momentum of the virtual (off-shell) emitter parton. We will also denote k = k 1 + k 2 , q = p − k. 4-dimensional transverse momenta k i⊥ = (0, k i , 0) for i = 1, 2 will be also used (to be extended to n = 4 + 2 dimensions wherever necessary). From k 2 1 = k 2 2 = 0 we obtain The lightcone variable decreases from 1 to x = 1−α 1 −α 2 after two emissions. The angular scale variable is a preferred choice, instead of transverse momentum. The virtuality of the emitter parton after two emissions (entering its propagator) and the gluon pair effective mass are: (3.4)

Inclusive evolution kernels, CFP vs. MC
In the CFP scheme, the NLO inclusive kernel is extracted from the second order expression for the bare PDF, which in compact CFP notation reads: where Z (2) F is the quark renormalization constant and K 0 = K [1] 0 + K [2] 0 is the 2-particle irreducible kernel (truncated to 2-nd order in perturbative expansion) defined in refs. [9,38]. The NLO contributions to the evolution kernels are coming from PK [2] 0 and PK [1] 0 ((1 − P)K [1] 0 ). Diagrams (b-g) in Fig. 1 are in the first class and only diagram (a) is in the second class (with subtraction). The other 2-nd order terms like (PK [1] 0 )(PK [1] 0 ) and Z [1] F PK [1] 0 yield pure 1 2 poles times elements of LO kernels and do not contribute to NLO kernels. The 1-st order 1Z [1] F + PK [1] 0 was already analyzed in the previous section. From now on we drop flavor indices as all but one diagram contributing to the non-singlet kernel at NLO level describes qq transitions. The qq flavor indices will be understood implicitly if not indicated otherwise. 12 The NLO contribution to the evolution kernel (to bare PDF of CFP) from a given 2R Feynman diagram X in Fig. 1 reads: where the θ-function limits phase space from the above using variable s(k 1 , k 2 ) = max{a 1 , a 2 } (resulting CFP kernel is independent of this cut-off), C is a color factor of a diagram X.
Monte Carlo featuring complete NLO evolution, can be expressed as a time-ordered exponential in the logarithm of its factorization scale Q, see [28]. Hence, at the inclusive level, it obeys its own evolution equation in ln Q with its own inclusive NLO evolution kernel, being the derivative in ln Q of the MC [28] distribution (truncated to 2-nd order) 13 Here K 0 is the same as in CFP and comes from the Feynman diagrams, albeit with realvirtual collinear cancellations executed before taking the derivative -so above formula is finally executed in n = 4.
For the use of P it is enough to apply eqs. (3.9,3.10) below. It acts on the integrand, contrary to P of CFP acting on the integrals, hence it provides unintegrated NLO distributions for the MC [20,28]. For our purpose (2R diagrams) the above reduces to the following: contains two-particle-irreducible diagrams and contains diagrams requiring soft counterterms. 14 The remaining action of P is spin projec- 12 Possible ways of implementing multi-flavor partons in MC have been presented in refs. [39,40]. 13 Strictly speaking, diagrams like (b-c) in Fig. 1 produce ln k (Q/µR), k > 1 (similarly to higher 1/ k poles in CFP), to be resummed into running coupling constant, before applying this formula. 14 Cut-off q0 plays no role in evolution kernel.
tion, the same way as in CFP. On the other hand, the subtraction term −P (K [1] 0 )P (K [1] 0 ) is identical to the double gluonstrahlung LO distribution of the MC, hence it deviates from CFP. For example interference diagrams contribute and subtracted NLO diagrams contribute (3.12) Since MC distribution W ct (k 1 )W ct (k 2 ) encapsulates (by construction) all collinear and soft singularities, subtracted G X a (Q, x) can be evaluated in n = 4. Alternative expressions for NLO inclusive kernels of eq. (3.6) and eq. (3.8) provide precisely the same results for graphs (d-e) in Fig. 1, which do not have internal divergences nor require subtraction (as in LO case), and are evaluated at n = 4. For diagrams (a-c) in Fig. 1 we shall see certain small but important differences, which indicate that the MC represents a slightly different factorization scheme than CFP.

Overview of the 2R phase space integration
It is convenient to introduce slightly differently normalized phase space 13) which is dimensionless in n = 4.
Most of the presented differential results will be normalized using the above integration element. For instance in eq. (3.6) we replace dΨ n (k) → dΦ n (k) and (3.14) Let us outline the general methodology used in the 2R phase space integrations. It will be described in n = 4, with some small modifications it will also apply in n = 4 + 2 . The integration procedure consists of the following steps: (a) Using the identity Θ Q>max{a 1 ,a 2 } ≡ Q 0 dQ δQ =max{a 1 ,a 2 } , the integration variableQ is introduced: 2π dφW X (a 1 /a 2 , φ, α 1 , α 2 ) Θ max{a 1 ,a 2 }>q 0 .
The first order expression for the PDF in the MC (performing scalar products) reads: (3.21) The above includes factors 1/2! due to Bose-Einstein (BE) symmetrization as well as 2 multiplying interference diagrams. Following the LO calculation example of Section 3.3, the integration over transverse degrees of freedom is done: Integration over α-variables finally provides: (3.23) As already said, P Bx (x) is the same in CFP and in MC schemes (up to a normalization factor 2), 16 because this diagram has no internal collinear divergence. Uncanceled IR divergences are still present (I 0 term). Summarizing, the distribution of eq. (3.19) will enter into the NLO correction to the MC exclusive kernel contribution. The distribution of eq. (3.23) will be used for numerical overall tests of the MC at the NLO level.

Subtracted double bremsstrahlung diagram
The double bremsstrahlung diagram of Fig. 1(a) (denoted as Br) is not 2PI (it consists of two 2PI LO kernels) and needs a subtraction term (referred to as diagram BrC).
We shall start with a simpler case of integrating Br−BrC contribution to evolution kernel in the MC scheme in n = 4. Next we shall recalculate the same Br−BrC contribution to the bare PDF and NLO kernel in the CFP scheme, analyzing all its components and discussing all differences with the MC case in detail.
The differential distributions for two Br diagrams 17 in n = 4 + 2 read as follows: where: (3.25) The most singular term ∼ T Br 2 can be rewritten (modulo O( 2 ) terms) as a product of two LO kernels: where z 1 = 1 − α 1 and z 2 = (1 − α 1 − α 2 )/(1 − α 1 ). The above term coincides in the MC for the ladder with the following counterterm, being just the LO MC distributioñ 27) It also enters as a subtraction term into the MC weight which implements NLO corrections. The contribution to the inclusive PDF of the NLO MC, including the explicit Bose-Einstein (BE) symmetrization factor 1/2!, reads: where is obtained from analytical phase space integration using the same methodology as in the CFP scheme, see below.
17 Two diagrams because of interchange of vertices due to the Bose-Einstein symmetrization.
It should be kept in mind that P Br sub (x) above is obtained for the angular ordering and it would be different if we would have adopted k T -ordering -the difference would be P Br Kin (x) of eq. (3.32) below, see discussion in refs. [33,37].
Let us now turn to the CFP scheme which provides, for this particular diagram, a 2R contribution to the NLO kernel, different than in the MC scheme. For calculating PK [1] 0 ((1 − P)K [1] 0 ) of the bare PDF of eq. (3.5) we follow procedure (a-f) of Section 3.3 step by step. After integrating overQ, φ and φ 2 , at step (e), we deal in the P(K [1] 0 K [1] 0 ) part with the singular integral in the variable y = max(y 1 , y 2 ): The most singular part due to the 1 δ y=0 term in P(K [1] 0 K [1] 0 ) is easily integrated: is just the double convolution of the LO kernel. The counterterm PK [1] 0 ((−P)K [1] 0 ) contributes the same expression, but with (−1) in front, and finally (PK [1] 0 )(PK [1] 0 ) adds the same expression with (+1) in front. Altogether, the pattern of building correct exponential structure of the LO in CFP is much more complicated than in the MC, with a lot of over-subtractions, see more discussion in ref. [29].
Our main aim however is the residue in front of the 1 pole. Most of it comes from the term 1 y + in eq. (3.30), which is in close correspondence with the NLO correction in the MC. In addition it gets "fall-out" contributions from 1 2 × terms. In particular x-independent terms from the expansion luckily cancel between P(K [1] 0 K [1] 0 ) and the collinear counterterm PK [1] 0 ((−P)K [1] 0 ). A similar but x-dependent contribution from T Br 2 gives rise to a remnant term The second bracket [. . . ] comes from the counterterm -it lacks 1 2! of the true distribution, and one of P (0) P (0) terms gets killed by the PP operation. Altogether, the above spin artifact of CFP (absent in MC) contributes the following: (3.31) The last important contribution in the class ∼ 1 4 δ y=0 × is due to the (α 1 α 2 ) 2 term in our particular choice of the phase space parametrization. In fact its role is to cancel the dependence on this choice, see also discussion in ref. [33]. It is produced by a similar mechanism of partial cancellation with the counterterm as above, and in the z-parametrization reads: .
Its contribution to the bare PDF and NLO kernel is: Finally the real physics is in the term 1 y + in eq. (3.30), which happens to be exactly the same (up to the normalization factor 2) as the MC contribution of eq. (3.28). At step (f) of the integration procedure it reads: (3.33) After α-integrations we obtain Γ Br sub = 1 2 P Br sub (x), (3.34) where P Br Sub (x) = 1 2 P Br Sub (x) of eq. (3.29). Altogether, the CFP kernel from the subtracted Br diagram is: reproducing the result of ref. [9]. As noted in ref. [33] P Br Kin (x) is absent in CFP, provided we choose the maximum transverse momentum as factorization scale variable: s(k 1 , k 2 ) = max(k T 1 , k T 2 ). It is simply the case because the factor (α 1 α 2 ) 2 is absent. However, an additional contribution exactly equal to P Br Kin (x) would appear from the 2R integral of eq.(3.28) due to adopting k Tordering, and it would exactly compensate the lack of P Br Kin (x) from (α 1 α 2 ) 2 . In a sense, the CFP scheme features an automatic self-correcting mechanism, such that it provides the result for k T -ordering, independently of the kinematic parametrization of the phase space actually used.
Summarizing on the subtracted Br diagram, the integrandW Br sub in eq. (3.28) will enter into the NLO correction to the MC distribution, while eqs. (3.29) and (3.35) will be used in the numerical tests of the MC codes. The difference between CFP and MC factorization schemes for the subtracted Br diagram (for Br+Bx as well) at the inclusive kernel is IR finite and reads: (3.36) Let us also sum up C 2 F contributions of the Br and Bx diagrams (using only C 2 F part of Bx) in the 2R phase space. For the MC scheme we obtain: (3.37) As we see IR part (I 0 term) cancels between Br and Bx diagrams. The same phenomenon is also true for the differential distributions, see figures and analytical investigation of the α i → 0 limit in ref. [32], which are based on the above results. The above IR cancellations are vital for the stability of the NLO MC weight, see tests of the prototype MC in ref. [28]. It should be stressed that it is not guaranteed 18 and here it is true thanks to a good choice of the multigluon LO MC distributions compatible with the soft (eikonal) limit already at the LO level, see eq. (3.27).

Gluon pair production diagram -Vg
Let us investigate now another important 2R NLO contribution from the gluon pair production diagram Vg of Fig. 1(b). We shall calculate its contribution to inclusive kernels focusing on possible differences between MC implementation and CFP scheme. As we shall see, the 2R gluon distribution from Vg diagram features very different singularity pattern than Br+Bx diagrams discussed previously -it has an internal collinear singularity for the mass of produced pair going to zero, a 2 → 0, that is located in a 1 → a 2 (instead of a 1 → 0).
Cancellation of this internal singularity happens without an intervention of (1 − P) operator, simply by adding a virtual diagram (gluon selfenergy). The remaining residual double logarithm (or double pole in ) singularity is related to the running coupling constant.
Mastering the additional soft-gluon sudakovian IR singularities α i → 0 will be again very important for stability of the NLO MC weight. A complete discussion of the soft gluon limit will require introducing the interference diagrams Yg and Bx of Fig. 1(d) and Fig. 1(e), hence it will be deferred until the next section.
From the Feynman diagram we obtain: (3.38) In the above we have kept only those terms O( 1 ) from the γ-trace which lead to 1 2 poles, because of the extra 1 from an internal 1 a 2 gluon mass singularity. Using the most singular term in eq. (3.38) is isolated even more clearly: (3.39) The term proportional to T V g 2+ ( ) = T V g + ( )+T V g 3 ( ) is the only one contributing the double pole 1 2 and is explicitly proportional to the product of two LO kernels: where z = α 1 /(α 1 + α 2 ) and P (0) . What enters into the 2R part of the NLO correction in the MC, see refs. [28,29], is not the above divergentW V g , but rather the non-divergent difference with the following "soft collinear counterterm" (SCC) representing the distribution used in the LO MC: where a max = max(a 1 , a 2 ). It is the result of a slight simplification of the term ∼ T V g 2+ in eq. (3.39). Generally, such a SCC is not unique, but complete NLO corrections to observables are insensitive to its choice. What is highly sensitive, however, is the dispersion (and positiveness!) of the MC weight implementing NLO correction. The above choice is unique in this sense, that it encapsulates not only collinear singularity from a 2 → 0, but also all soft gluon singularities α i → 0 for all three diagrams Vg+Yg+Bx, see next section.
This property ensures good behavior of the NLO MC weight. Moreover, the factor in eq. (3.41) can be iterated into a separate final state LO sub-ladder for the gluon emitted from the primary initial state ladder, see refs. [28,29].
As we are working on the exclusive level we are technically similar to the techniques of hard process subtractions of ref. [41] (dipole subtraction) or ref. [42] (antenna subtraction), 19 but not to the subtractions used in the inclusive calculations of NNLO evolution kernels of refs. [10,11].
In view of the above discussion, it is useful to know analytically, for numerical crosscheck of the MC code and for discussing complete NLO corrections in the ladder MC, the following subtracted Vg contribution to the NLO PDF 20 calculated in n = 4 (as usually including BE factor) (3.42) Generally, the presence of the SCC subtraction is a natural element in any MC scheme with soft gluon (photon) resummation, either to the hard process or to the ladders, in order to eliminate possible double counting of the singular term. However, the use of subtractions can also simplify non-MC calculations, like analytical integration of the Vg diagram over the 2R phase space in n = 4 + 2 dimensions in the CFP scheme, before combining it with the gluon self-energy virtual diagram. Let us comment more on that, venturing a little bit in the area of combining 2R and 1R1V contributions (complete discussion is beyond the scope of this work). In such a case, it is useful to split 2R gluon phase space into a > κa max and a < κa max , κ 1, schematically and split the Vg contribution into a subtracted one and the SCC. The subtracted part in the decomposition 19 The dipole/antenna subtractions works for at least two ladders (they are dealing with hard process), whereas our method works within one ladder. Furthermore, we pay special attention to the fact that our counterterms can be iterated in the MC simulation. 20 In the above 1/αi are regularized by a small parameter δ as in CFP, however, this is not necessary once diagrams Yg and Bx are added, see next section.
can be evaluated in n = 4, all over the phase space. On the other hand, the SCC part Γ CT is evaluated analytically separately in the "resolved" part a > κa max in n = 4, and separately in the a < κa max part in n = 4 + 2 . This allows to profit from adjusting phase space parametrization to specific complications of the integrand in each part! Finally, one combines all three parts into a formula for the 2R integrated Vg in n = 4 + 2 . The parameter κ, and even the dependence on the particular choice of SCC drops out in the final result. The other immediate profit from the above methodology is that one may combine 2R from a < κa max with the 1R1V contribution (gluon self-energy) such that the Sudakovian part of the 1 2 pole gets cancelled. 21 This kind of calculation for Vg in n = 4 + 2 is presented in the Appendix using a slightly different choice (for historical reasons) of the SCC: (3.43) Switching from one kind of SCC in the integration to another is relatively simple, see Appendix.
Last but not least, let us discuss the differences between the MC scheme and the CFP scheme for the Vg diagram. In the previous case of Br diagram we have seen that the basic mechanism of producing differences between two schemes is the action of the 1 − P operator. Since this operator is absent for Vg, one generally expects no differences between the two schemes. In particular any effect of the terms proportional to from γ-traces will land in the ∼ δ(a 2 ) part, where soft 2R and 1R1V are combined together in the same way in both schemes.
The only subtle point is the term b 0 1 2 left over from adding 2R soft and 1R1V contributions. It comes from integrating over the ln q µ R term in the gluon self-energy, and in the MC scheme it builds up an α s dependence for some kinematic variable. Which variable? Changing from one choice of q to another may generate an extra NLO term in the kernel (in MC scheme). Our preliminary study shows that taking transverse momentum as q is compatible with CFP, that is the Vg diagram contribution is then identical in MC and CFP. The contribution from the diagrams Yg and Bx discussed in the next section, will contribute the same way in both schemes, due to the lack of any internal collinear singularities.

Gluon interference diagram -Yg
The Monte Carlo distribution for the gluon interference diagram of Fig. 1(d)  where: The expression for the contribution of the Yg diagram to the NLO kernel is equal to (we include both factors 1/2! from BE and 2 due to interference): The integrand only has the familiar scale singularity which can be extraced in standard way as a logarithm. Then the integrals over transverse components take form: .

(3.47)
After the final integration over α i variables we obtain , (3.48) Gluon interference diagram Yg features both double and single logarithmic IR divergences (I 1 and I 0 ). Of course, they must cancel when all diagrams are added. Cancellations occur for 2R diagrams not only on the inclusive integrated level but already on the exclusive unintegrated level, see [32]. We can see them explicitly by adding inclusive kernel contributions for gluon interference diagram Yg (3.48), subtracted gluon pair production diagram Vg (3.42) and part of gluonstrahlung interference diagram Bx (3.23) proportional to 1 2 C A C F colour factor: (3.49) The above expression for the sum of 2R contribution with colour factor equal 1 2 C A C F is free from double and single logarithmic IR divergences (I 1 and I 0 ), as expected.

Other non-singlet diagrams
The contributions to NLO kernels from diagrams displayed in Fig. 1(c), Fig. 1(f) and Fig. 1(g) will be presented in the following. They will be referred to as Vf, Yf and Xf diagrams respectively. As discused in the introduction, before all singlet class diagrams are included, the contributions from the interference diagrams Xf and Yf should enter the MC code in the inclusive form.

Interference diagram Xf
The crossed-ladder diagram Xf contributes the following distributioñ to the quark-antiquark kernel, where: , T Xf 12 = 2(x 2 + 1) The integrated distribution is equal to: and the contribution to the NLO kernel equals: The above agrees with [9] up to the − sign which is a matter of convention.

Fermion interference diagram -Yf
The differential distribution for the fermion interference diagram Yf of Fig. 1(f) reads: where: (4.6) The contribution of the Yf interference diagram to the inclusive kernel (bare PDF) is given by (including interference factor 2): , ) .

(4.7)
This result agrees with that of ref. [43], 22 the difference in sign can be attributed to a different definition of space dimension (n = 4 − 2 , being negative as opposed to CFP).

Fermion pair production diagram -Vf
The fermion pair production diagram Vf of Fig. 1(c) features an internal collinear singularity when the mass of the produced pair goes to zero. The kinematical structure of the Vf graph is quite similar to that of the gluon pair production diagram Vg. The n = 4 + 2 dimensional distribution for the Vf diagram reads: (4.8) 22 The authors use a different regularization technique. In case of the Yf diagram, however, both methods give the same result, as explained in [43].
-25 - where: (4.9) The contribution of this diagram integrated in n = 4 + 2 dimensions reads: (4.10) For exclusive modeling of this diagram we need to deal with its internal collinear singularity in a similar way as for the gluon pair production diagram Vg. The following counterterm can be used both for the MC purpose in n = 4 and for combining real and virtual contributions. In the latter case, we decompose the Vf contribution into a part entering Monte Carlo a>κamax and an unresolved part Γ V f CT a<κamax required for the cancellation of double poles from the virtual contributions. For completeness we give the subtracted contribution of the Vf diagram to the NLO PDF: , (4.12)

Summary and outlook
The main result of this work is a complete collection of 2-real parton (quark, gluon) differential (unintegrated) distributions, which enter calculations of the NLO DGLAP non-singlet evolution kernels, in a form ready for the use in the Monte Carlo implementation of the ladder (also referred to as NLO parton shower MC). The distributions are given in eqs. (3.19), (3.24), (3.38), (3.44), (4.1), (4.5), (4.8). These distributions in the fully differential form are not available in the literature. We also present the differential distribution of the collinear soft counterterms, which are used to subtract internal singularities for some diagrams. These subtractions are also present in the MC weights. The MC collinear soft counterterm distributions are defined in eqs. (3.27), (3.41).
Furthermore, we present analytical integration results. They are presented as the contributions to evolution kernels from the same 2-real parton differential distributions listed above, see for example eqs. (3.37) for all bremsstrahlung diagrams. In case of diagrams with internal collinear singularities, subtractions of the MC collinear counterterms is done. For certain diagrams it was possible to compare the integration with available published results of refs. [9,43,44] and agreement was found.
The QCD evolution of the NLO ladder implemented in the MC is slightly different from that of standard M S, as defined and implemented in Curci-Furmanski-Petronzio paper [9]. For instace, the differences between MC and CFP schemes are discussed as far as it is possible for the 2-real contributions. They are typically present in the diagrams with internal, collinear divergences, see for instance eq. (3.36). The complete discussion of this issue is beyond the scope of the present work -it will be completed when diagrams with 1 real and 1 virtual corrections are added into the game in the forthcoming work. Nevertheless, even incomplete results provide us important insight into the differences between NLO (integrated) kernels of MC and CFP M S schemes. 23 This analysis will also be a practical outcome of the entire project.
We did not explicitly show results of the numerical cross-checks of the analytical results. Let us only mention that all analytical integration results in the paper were cross-checked (up to 4-digits) by means of the MC numerical integration using FOAM program [45] within the MCdevelop system [46].
Summarizing, the present work marks an important step forward on the way to the implementation of the complete NLO DGLAP ladder in the Monte Carlo form.

A.2 Two counterterms for Vg
For various purposes it is useful to calculate the contribution from both counterterms of eqs. (3.43) and (3.41) in n = 4 with the cutoff |a| > κa max (ã > κ). Let us start with the 24 We always implicitly use Principal Value type of regularization for α integrals.  Figure 3. The disc represents the region |a 1 | < |a 2 |. The dimensionless ratio of y = |a 1 |/|a 2 | and φ, the relative angle between a 1 and a 2 , are shown. Alternative variables (ã, θ) are also shown.
counterterm of eq. (3.43). Introducing dimensionless variables y i = a i /a max ,ã = |a|/a max and performing integration over the overall scale a max we obtain: We need to remember that the κ cutoff is infinitesimal and all terms O(κ) are neglected. Both y 1 and y 2 integrals are equal becauseã 2 (y 1 , y 2 ) is symmetric andq 2 (1, 1) = (1−x) 2 α 1 α 2 does not depend on y i .
α S 2π 2 dα 1 dα 2 α 1 α 2 δ 1−x=α 1 +α 2 × 2π 0 dφ 2π dy 1 y 1 T V g + a 2 (y 1 , 1) + T V g 3 (y 2 1 − 1) 2 a 4 (y 1 , 1) θã >κ . (A.5) The integration over α factorizes from the y and φ integrals, hence it is performed separately. Moreover, because of the cutoff onã 2 it is convenient to change variables. Instead of (y, φ) variables we use (ã, θ) depicted in Fig. 3. The jacobian for this transformation is equal toã/ √ 1 +ã 2 − 2ã cos θ. The two integrals above are calculated separately. Firstly, 1 2π Combining both parts together we obtain: (A.8) After performing the α-integration the final result for the first counterterm of eq. (3.43) reads: (A.9) The integration for the counterterm of eq. (3.41) representing the double gluon distribution in the Monte Carlo proceeds for the same phase space quite similarly and the final results reads: The difference between the two counterterms is, of course, colliner-convergent and it reads P K−CT (x) = P K a>κ (x) − P CT a>κ (x) = α S 2π The above is used to correct the Vg subtracted result of eq. (A.3) in order to obtain the result of eq. (3.42).