Automating dipole subtraction for QCD NLO calculations

In this publication the construction of an automatic algorithm to subtract infrared divergences in real QCD corrections through the Catani-Seymour dipole subtraction method [arXiv:hep-ph/9605323] is reported. The resulting computer code has been implemented in the matrix element generator AMEGIC++ [arXiv:hep-ph/0109036]. This will allow for the automatic generation of dipole subtraction terms and their integrals over the one-parton emission phase space for any given process. If the virtual matrix element is provided as well, this then directly leads to an NLO QCD parton level event generator.


Introduction
Perturbative calculations form one of the best understood methods to provide predictions for the behavior of a Quantum Field Theory and to compare them with experimental results. Many of the methods applied in such calculations have found their way into textbooks already decades ago, see e.g. [3]- [7]. Typically, the perturbation parameter is related to the coupling constant of the theory in question, which in most cases indeed is a small quantity. This also implies that the corresponding fields may asymptotically appear as free fields and thus are the relevant objects of perturbation theory. Obviously, this is not true for the strong interactions, i.e. QCD, where the fields, quarks and gluons, asymptotically are confined in bound states only. This is due to the scaling behavior of the coupling constant of QCD, α S , which becomes small only for large momentum transfers, see for instance [8,9]. It is the confinement property that to some extent restricts the validity of perturbative calculations in QCD to the realm of processes characterized by large momentum transfers or by other large scales dominating the process, such that factorization theorems can be applied [9]- [11].
Typically, for most of the relevant observables in particle phenomenology, the leading term of the perturbative expansion can be related to tree-level diagrams. In the past years, the calculation of these terms has been fully automated and a number of tools capable of dealing with up to eight to ten external particles without any significant user interference have emerged [2,12]- [16]. However, for many practical purposes, tree-level calculations are not sufficient. This is due to a number of reasons: first of all, many measurements aim at the extraction of fundamental parameters. However, in Quantum Field Theories, parameters are subject to corrections, which usually exhibit ultraviolet divergences. These divergences are dealt with through the renormalization procedure, which can be done in a scheme-and scale-dependent way only, see e.g. [6,17]. Therefore, in order to extract parameters from the comparison of a (perturbative) calculation with experimental data, the calculation itself must contain the same kind of quantum corrections necessitating their renormalization. Second, it should be stressed that in tree-level calculations, there are some choices to be made, concerning the scale at which inputs such as the coupling constant, quark masses or parton distribution functions are taken. In principle, different scale choices are equivalent, and renormalization group theory guarantees that, when taking into account all orders, the effect of scale choices vanishes. At leading order (LO), however, their impact may still be significant, such that tree-level calculations merely give the order of magnitude for corresponding cross sections etc.; a prime example for this is the production of a Higgs boson in gluon fusion processes, where only the next-to-next-to leading order correction significantly reduces the scale dependence and produces a stable result [18,19]. Thus, aiming at any more precise prediction, higher-order calculations are a crucial ingredient of phenomenological analyzes. But although indispensable, so far there is no fully automated tool available for QCD calculations at next-to leading order (NLO), i.e. at the one-loop level. This is because a true NLO calculation is certainly much more complex than a leading order (LO) one. First of all, some of the essential ingredients, namely the loop or virtual contributions are not under full control yet. In general, up to now calculations of these corrections to physical processes are limited to contributions containing five-and in some cases six-point functions, see for example [20]- [25]. But even to reach the level of known scalar master integrals is far from being trivial; the tensor reduction necessary for this step [26] results in a proliferation of terms with non-trivial cancellations among them, which render the implementation in a computer code a major effort. On the other hand, some of the loop corrections exhibit not only ultraviolet divergences to be renormalized, but also infrared divergences. They also need to be regularized, but then they must be canceled against similar infrared divergences stemming from the real contributions. This basically translates into canceling divergences in phase space volumes of different dimensionality. The cancellation in fact is one of the most important consequences of the Kinoshita-Lee-Nauenberg or mass factorization theorems [27]- [29]. However, in order to practically achieve the cancellation, the real infrared divergences also need to be regularized. Essentially, there are two ways of doing this.
One method, also known as phase-space slicing [30]- [35], bases on dividing the phase space of the additional real emission into an infrared-safe (hard) and a infrared-divergent (soft) region. The division is usually performed by subjecting pairs of particles to an invariant mass criterion. Then, the soft region is integrated analytically in d dimensions. Typically, in this step, the helicitysummed matrix element squared is approximated by its double-pole (or eikonal) limit. The result of the analytical integration will contain single or double poles of the form 1/(d − 4) or 1/(d − 4) 2 , respectively. They typically are accompanied with logarithms of the invariant mass criterion. Such logarithms, but with opposite sign, also appear in the numerical evaluation of the full matrix element squared for real emission in the hard region of phase space, performed in 4 dimensions. In principle, these two potentially large contributions (logarithms of a potentially small quantity) originate from the unphysical division of the phase space and should thus cancel. Therefore, the key issue thus in phase space slicing is to adjust the parameters of the procedure such that the dependence on the slicing parameter is minimized. So far, this adjustment has been done manually only and this is one of the reasons why other methods have become more popular with practitioners of NLO calculations. Such alternative methods of dealing with the real infrared divergences base on directly subtracting them [1,36]- [42] 1 . At NLO level the subtraction in all methods is performed such that the additional particle is added to the leading order matrix element in a well-defined way through terms which, on one hand, exhibit the correct divergent behavior in the soft and collinear limit, and, on the other hand, can easily be integrated over the full d-dimensional phase space of the extra particle. The idea is then that the so subtracted matrix element squared is finite and thus can safely be integrated numerically in 4 dimensions. On the other hand, the subtraction term is added to the virtual bit and integrated analytically in d dimensions. Again, it exhibits single or double poles of the form 1/(d − 4) or 1/(d − 4) 2 , respectively. These poles again cancel the infrared poles of the virtual contributions. The fact that there are universal subtraction terms, i.e. terms which will cancel the infrared divergences in a process-independent manner, is one of the main reasons why subtraction methods have become increasingly popular in past years and why they have been used for many of the state-of-the-art calculation of NLO corrections to physical processes, like for instance [20]- [24].
The universality of the subtraction terms also allow for an automated treatment of real infrared divergences. It is the subject of this publication to report on a fully automated, process-independent implementation of one of the popular subtraction procedures, ready for use in realistic NLO calculations. Therefore, the outline of this paper is as follows: In section 2, the anatomy of QCD NLO calculations will be formalized in a more mathematical language and the chosen subtraction method, the Catani-Seymour dipole subtraction [1] will briefly be reviewed in its original form for massless particles. Although its extension to massive particles [40] is straightforward from an algorithmic point of view, this paper concentrates on the massless case only. In section 3, the fully automated implementation of the corresponding massless dipole subtraction of arbitrary matrix elements into the matrix element generator AMEGIC++ [2] will be presented in some detail. Some simple tests of the implementation will be discussed in section 4, before some physical applications and the comparison with results from the literature will round of the presentation in section 5.
2 Brief review of the Catani-Seymour formalism 2.1 NLO cross sections and the subtraction procedure Cross sections at NLO precision are given by where the LO part σ LO is obtained by integrating the exclusive cross section in Born approximation over the available phase space of the m final state particles and, eventually, over the Bjorken-x of incident partons. Ignoring this additional complication for the sake of a compact notation, The LO cross section is thus given by where d (4) Here, d (4) Φ (m) denotes the phase space element of m particles, taken in four dimensions, M m is the matrix element for the process under consideration, and F (m) J is a function of cuts defining the jets etc.. As already indicated, here and in the following, the superscripts in the integral denote the dimensionality of the integration. In order to obtain a meaningful result to be compared with experimental data, typically isolation cuts are applied on the outgoing particles, which may also serve the purpose of keeping the integral finite. A typical criterion for example is to identify outgoing partons with jets and thus apply jet definition cuts on the partons such that they are all well separated in phase space. Anyway, the cuts will not be stated explicitly in the integral, but they are understood implicitly with the integration, including suitable generalisations in d dimensions, where necessary. Thus, the integration of the Born level cross section can directly be carried out in four space-time dimensions, as indicated in the equation. In view of the dipole subtraction formulae, it is useful to introduce at this point bras and kets m 1, . . . , m ′ | and |1, . . . , m ′ m . They denote states of m final state partons partons labeled by 1 to m ′ and are vectors in colour and helicity space. Introducing, in a similar fashion, vectors for the spins and colours, matrix elements thus can be written as Therefore, in this notation, the matrix element squared, summed over final state colours and spins reads The two integrals on the right-hand side of Eq. (6) are separately infrared divergent in four dimensions, and are therefore taken in d dimensions. For the real correction, the divergences arise when the additional parton becomes soft or collinear w.r.t. some other parton, leading to on-shell propagators in the matrix element. For the virtual correction, the divergence comes with the integration over the unrestricted loop momentum, such that again a propagator goes on-shell. As already stated in the introduction, now the celebrated theorem of Kinoshita, Lee and Nauenberg [28,29] comes to help and guarantees an exact cancellation of two divergent contributions, thus keeping their sum finite 2 . Setting d = 4 + 2ǫ in the following, the divergences will manifest themselves in double and single poles, i.e. as 1/ǫ 2 and 1/ǫ, respectively. In principle, cancellation of the poles then solves the problem; in practice, however, the direct applicability of the equations above to real physical processes is limited since analytical integration over a multi-particle phase space in d dimensions with cuts in many cases is beyond current abilities. Therefore, a detour has to be taken. The idea is to construct a subtraction term for the real emission contribution, which encodes all of its infrared divergences, but can analytically be integrated over in d dimensions. In this way the infrared pole structure of the real part with its 1/ǫ and 1/ǫ 2 poles is exhibited and cancels the corresponding virtual contributions. Subtracting this term from the real emission contribution and adding it to the virtual corrections then eliminates the infrared divergences in both parts. The subtracted real matrix element squared then is finite and thus its full (m + 1)-particle phase space can safely be integrated over in four dimensions. In this way, the subtraction term aims at an infrared regularisation of the two contributions at integrand level.
The catch of the subtraction method now is that the subtraction terms can be obtained from the Born terms in a straightforward way and that only the phase space integral of the extra particle has to be taken in d dimensions, while the phase space for the remaining m particles can be taken in four dimensions. This is similar to the way, the loop terms are evaluated. There, only the loop integration is performed in d dimensions, whereas the phase space of the outgoing particles is done in four dimensions. Therefore, the final structure reads Both integrands now are finite, allowing all integrations to be performed numerically. In contrast to some other regularisation methods (like, e.g., phase space slicing) the subtraction method does not rely on any approximation and does not introduces any ambiguous and/or unphysical cut-off scales etc., as long as the integration of d (d) σ A can exactly and analytically be performed.
In [1] a general expression for d (d) σ A has been presented, called the dipole factorisation formula, allowing to write 2 In fact, this is only guaranteed for infrared-safe quantities. More specifically, if F (n) J defines jets in terms of the momenta of an m-parton final state (taken at Born level), infrared safety demands that F such that, symbolically, where Here the sum of the dipole terms V dipole contains all soft and collinear divergences of the real emission pattern. This factorisation formula is suited for any process with massless partons, and fulfills all the requirements mentioned above. An extension to massive partons has been presented in [40].
However, as already mentioned in the introduction, in this publication only the massless case will be considered. In order to provide a self-contained description, all necessary analytic expressions will be listed in this publication.

Generalisation to hadronic initial states
The cross sections discussed so far were given for point-like initial states. For cross sections in hadron collisions, however, the differential cross sections above must be convoluted with parton distribution functions (PDFs): Here the subscripts on the cross section denote the flavours of the incoming partons; for the total cross section a sum over them has to be performed. For the NLO part, now the higher-order corrections residing in the PDFs must be taken care of. This is done by supplementing the NLO part with a collinear subtraction term dσ C ab , such that This new term contains collinear singularities, incorporated in 1/ǫ-terms and reads The collinear subtraction term is factorisation-scale and scheme dependent. This scheme dependence resides in the terms K F.S. , which, for the common M S-scheme vanish, i.e. in this scheme all terms K F.S. = 0. However, this scheme dependence cancels similar terms in the PDFs such that, taken together, the full hadronic cross section again is scheme-independent.
In the case of incoming hadrons, the subtraction method is applied to σ NLO (p a , p b , µ 2 F ) as described before, with the only difference that in this case the singularities of dσ V ab only cancel in the sum

Observable-independent formulation of the subtraction method
Up to now, the dσ denoted cross sections in a broad sense. To be a bit more specific consider the following expression for a cross section at Born-level and the corresponding next-to leading order expression: where dΦ (n) represents an n-particle phase space element, and M (m) , M (m+1) and V (m) are the LO matrix element, the NLO real matrix element and the NLO virtual correction matrix element, respectively. F (n) is a function that defines a cross section or an observable in terms of the nparton momentum configuration. In general, the function F may contain θ-functions (to define cuts and corresponding total cross sections), δ-functions (defining differential cross sections), kinematic factors or any combination of these. However arbitrary this sounds, there is a formal requirement on this function F , namely that in the soft and collinear limits, i.e. for cases where one parton becomes collinear w.r.t. another one or where one parton becomes soft, the function F (m+1) reduces to F (m) : The first two conditions define infrared-safe observables -to phrase it intuitively this means that such infrared-safe quantities must not be altered by additional soft or collinear activity. The last condition above is required to properly define the Born cross section.
Applying the subtraction method to the NLO-part of Eq.(16) results in where d[p] is the phase space element for the 1-parton phase space. In order to have an identity between the subtracted terms and the added term, both the (m + 1)parton contribution and the m-parton contribution have to be subjected to the same function F . To be able to perform the integration over the one-parton phase space independent of the observable this function therefore must be F (m) . In the case of the (m + 1)-parton contribution F (m) is applied to the m-parton configuration, generated by corresponding mapping given in the prescription of the dipole function.

The dipole subtraction functions
The universality of the soft and collinear limits of QCD matrix elements are the basis for the construction of the dipole subtraction terms. In both limits any matrix element squared for m + 1partons factorizes into an m-parton matrix element times a (singular) factor.
To be specific, consider first the soft limit of the matrix element, given by the momentum p j of parton j becoming soft, i.e. p µ j = λq µ with λ → 0 []. Then, employing the soft limit reads m+1 1, . . . , j, . . . , m + 1|1, . . . , j, . . . , m + 1 m+1 In a similar way the limit where two partons i and j become collinear is defined through p j → (1 − z)/z p i . In this limit the (m + 1) parton matrix element can be rewritten as where, again, theP (ij),i (z, k ⊥ ) are the well-known Altarelli-Parisi splitting functions. Then, the actual dipole function generating the limit, where one of the partons i, j of a m + 1-parton configuration becomes soft or both partons become collinear to each other, symbolically has the following structure: with the non-singular m-parton matrix element m ...||... m and the operator V ij,k , describing the splitting of the parton (ij). Here, and in the following, the splitting kernels V ij,k are matrices in the helicity space of the emitter. The dipole function also involves a third parton as 'spectator'. This parton in fact is identical with the colour partner k in the soft limit, Eq. (20). The form of the subtraction means that kinematically, 3 → 2 mappings are considered such that all involved partons are allowed to remain on their mass shells.
In general the splitting parton (called 'emitter') and the spectator can be both, initial and final state particles. This discriminates four different types of dipole functions, displayed in Fig. 1. The full subtraction term for any matrix element with (m + 1) partons in the final state is given by the sum of all possible dipole functions. For the most general case with two partons in the initial state, therefore In the following the explicit expressions for the dipole functions will be listed. The corresponding one-parton phase space integrated subtraction terms are discussed in Sec. 2.5.

Final state emitters with final state spectators
The dipole contribution D ij,k for the singular limit p i ·p j → 0, where all three involved partons are in the final state, is given by It is obtained from an (m + 1)-parton matrix element by replacing the partons i and j with a single partonĩj, the emitter, and the parton k is replaced byk, the spectator. The flavours of emitter and spectator are assigned as follows: The spectatork remains unchanged, and the emitterĩj is defined by the splitting processĩj → i + j. The product of colour charges in the numerator of Eq. (25) introduces an extra colour correlation in the m-parton matrix element. The kinematics of the splitting are described by the following variables and to obtain the momentaĩj andk in the m-parton configuration the following map is being used: Obviously, four-momentum conservation is exactly fulfilled, i.e.
and all partons remain on their mass shell, The splitting matrices, which are related to the d-dimensional Altarelli-Parisi splitting functions, depend on the spin indices of the emitter parton. For the case of a quark splitting (using helicity indices s and s ′ ) the kernel is a matrix in helicity space, whereas for gluon splittings (to a quarkanti-quark pair or to gluons), the splitting matrices are given by Lorentz tensors. This yields respectively. The dipole terms given in this section are sufficient for the subtraction procedure in the case of non-hadronic initial states such as e − e + -annihilation.

Final state emitters with initial state spectators
For the case of an emitting final state parton, the presence of an initial state spectator results in additional contributions to the singular limit p i ·p j → 0 of the full m + 1-parton matrix element. The corresponding dipole terms in this case are given by The kinematic variables now read and the momenta of the m-parton configuration are obtained by the map Again, four-momentum conservation is trivially fulfilled and the partons remain massless.
The corresponding splitting functions used in Eq. (31) read

Initial state emitters with final state spectators
The next type of dipole function now covers initial state singularities p a · p i → 0 with final state spectators, given by The partonãi, which enters into the m-parton matrix element on the r.h.s. of Eq. (35) is given by the splitting of the initial state parton a →ãi + i. The relevant kinematic variables in this case are and the momenta for the m-parton configuration are obtained bỹ The splitting matrices V ai k in Eq. (35) are The three dipole types discussed up to now (FF, IF, FI) are sufficient to construct the subtraction term dσ A for processes with exactly one initial state parton, i.e. DIS configurations.

Initial state emitters with initial state spectators
The remaining dipole function, only required by processes with two initial state partons, covers the case where both, the emitter and the spectator, are initial state particles, To describe the splitting, the following kinematic variables are used The construction of the m-parton kinematics for this dipoles differs from the other three cases. The reason is that in this case the emitter and the spectator are fixed to remain along the beam axis. Therefore all final state momenta (not only momenta of QCD partons) are transformed according to the map where The momentum of the spectator p b remains unchanged. The transformation above can also be interpreted as applying a rotation and a boost turning initial state momenta back to the beam axis after a mapping similar to the first three cases of dipole functions. Indeed it can be shown that the transformation of final state momenta in Eq. (41) is just a Lorentz transformation. However, in this case, the splitting matrices read

Phase space factorisation
In order to combine the poles of the subtraction function and the virtual matrix element the subtraction function has to be integrated analytically over the one-parton phase space of the respective splitting. The rules for the momentum mapping from 3 to 2 parton phase spaces have been constructed in Secs. 2.4.1-2.4.4 such that the corresponding phase space exactly factorizes.
As an example, and in order to fix the notation, the case of a final-final dipole, D ij,k , will be discussed in the following. There, the three-particle phase space for the partons i, j and k (all other partons are not affected by the splitting and will be omitted) in d dimensions is given by This can be factorized in terms of the mapped momenta, such that where [dp i (p ij ,p k )], written in terms of the kinematic variables defined in section 2.4.1, reads Within the dipole function only the splitting function itself depends on the variablesz i and y ij,k . Thus, the integration in d dimensions can be performed once and for all, independent of the specific scattering process under consideration. The result of the integration for each splitting type can be expanded as a Laurent series including double poles (∼ 1/ǫ 2 ), single poles (∼ 1/ǫ), and finite terms (∼ ǫ 0 ). Further terms of O(ǫ) are unimportant here and will be left out.
All results for the final-final and for all other dipole types can be found in [1].

Full result
Having at hand the integrals for each dipole function, all individual dipoles present in a specific process can be collected to yield the overall infrared divergence of the subtraction term. Then, the starting point for the calculation of jet cross sections in the dipole subtraction formalism reads where {m+1} denotes the sum over all parton-level processes. However, the important point here is to exactly cancel the poles of the corresponding individual one-loop parton-level processes, which is done exclusively for each momentum and flavour constellation. Therefore, for each specific m-parton process at NLO only a selection of dipole functions related to (m + 1)-parton processes contributes to the cancellation of the virtual divergences. In [1] it has been shown that this amounts to an effective reordering of phase space integrals and sums over parton configurations, such that where dσÃ {m} is the integrated dipole term that collects the integrals of all dipole functions and thus cancels the singularities of dσ V {m} . It is explicitly given by where dσ B {m} × I(ǫ) is a shorthand for the following procedure: Write down the expression for dσ B {m} , and replace the corresponding squared Born-level matrix element using the insertion operator I(ǫ) as defined below.
Finally, the full result for the integrated dipole term and the collinear counterterm as defined in Eq. (14) for the most general case with hadronic initial states reads where a and b again specify the initial state partons. The summation over a ′ and b ′ runs over all parton flavours, i.e. it includes gluons, quarks and anti-quarks occurring in the PDF.
The insertion operator I reads where the indices I and J run over initial and final state partons. The universal singular functions V I (ǫ) depend merely on the flavour of I and are given by with N f being the number of contributing quark flavours. The complete singular structure in Eq. (52) is contained in [dσ B ab (p a , p b )×I(ǫ)] and the sum [dσ B ab (p a , p b )× I(ǫ)] + dσ V ab (p a , p b ) must be finite for ǫ → 0. The finite insertion operators K and P are given by and Note that here the index i runs over final state partons only. The flavour-dependent functions K aa ′ (x),K a,a ′ (x), and P a,a ′ (x) are defined in Appendix A. As already mentioned, the factorisationscheme dependent function K aa ′ F.S. (x) vanishes in the commonly used MS-scheme. To obtain the final result for processes with no initial state partons only the I-term needs to be considered in Eq. (52). For processes with one initial state parton only, the result is obtained by using the I-term and one of the two integrals over K and P only, while omitting the contribution of K a,a ′ (x).

Freedom in the definition of dipole terms
As stressed before, the singular limits of the dipole functions are fixed by the requirement to cancel the singularities of the real correction matrix element. However, away from this limit there is some freedom for modifications. One possible modification has been presented in [21], where a parameter α has been introduced which cuts off a dipole function for phase space regions far enough away from the corresponding singularity. Its main advantage lies in a significant reduction of the average number of dipoles terms to be calculated for each phase space point of the (m + 1)-parton phase space of the real correction term. This constitutes an important alleviation of the calculational burden, since the total number of dipole terms grows approximately as m 3 . The α-modified subtraction terms also allow nontrivial checks of the implementation, since the total result must be independent of α. The α-modified dipole functions have been defined as follows: They will be employed later, in the implementation presented in this paper. Of course, such a redefinition of the splitting kernels also requires a recalculation of their integrals. The new αdependent insertion operators I and K have been presented in [21].
Another simple modification is the addition of finite terms to the splitting functions, such as The constant C directly ends up as a finite term in the integral of the splitting function and thus it can be easily included in the insertion operators of I and K, too. This again allows checks of the implementation, but it can also be employed to improve the numerical behaviour of the phase space integrals and to reduce the number of negative events.

Implementation in AMEGIC++
The Catani-Seymour dipole subtraction terms have been implemented in full generality into the automatic matrix element generator AMEGIC++, based on it's version 2.0 [52] 3 . In particular this translates into AMEGIC++ being able to automatically generate all relevant parts of the NLO matrix element within the subtraction method except for the virtual matrix element. It can be applied to any process with massless partons for which the real correction ME can be generated, an extension to allow also for massive particles is foreseen. This includes standard model processes as well as implemented extensions, as long as there are no new strongly interacting particles involved. For standard model processes the boundary is currently at about six-eight partons (initial and final state). The new implementation aimed at a maximal reuse of already developed automated methods of amplitude generation and process management. Therefore, first a brief overview over the relevant parts of the code are given before the new implementation is described in some detail.

Amplitudes
The matrix element evaluation in the C + +-code AMEGIC++ is based on the evaluation of Feynman amplitudes using a helicity method based on the developments in [57]- [59]. The fundamental idea of this method is to introduce a helicity basis, in terms of which all physical spinors can be expressed. This allows to compute each amplitude directly as a complex function of physical momenta and helicity/spin states instead of computing traces of spinor products and γ-matrices for squared amplitudes. The colour within any amplitude is treated separately, i.e. in the first step all colour factors (the SU(3) structure constants f abc and t a ij ) corresponding to the k-th amplitude A k are collected in an array C k .
Squared matrix elements can thus be written as and hence a colour matrix of complex numbers is generated once for each process. In an initialisation run, AMEGIC++ generates all formulae for the amplitude calculation of the user-specified parton level processes including the colour matrix c i,j (using a set of replacement rules for the colour algebra), it simplifies the expressions for the helicity amplitudes by identifying and factoring out common pieces and finally stores everything in C++ libraries.

Process structure and organisation
Since typically many parton level processes contribute to jet cross section calculations, usually a long list of processes needs to be computed. The corresponding structure in AMEGIC++ is as follows: • Any parton level process is represented by the class Single Process, • Process Group contains a (possibly recursive) list of such processes or groups of processes.
All parton level processes sharing a specific common set of properties are grouped together in two or three levels of groups. In many cases there are subprocesses contributing to the same jet cross section which are very similar. Therefore AMEGIC++ applies a procedure to identify such processes in order to save computer resources and accelerate the calculation. The following checks are performed: • Direct comparison of amplitudes: check for processes that have identical graphs, where all involved particles have the same masses, widths and underly the same interactions (with coupling constants that differ at most in a constant factor). Example: QCD processes that differ in quark flavours only.
• Numerical comparisons: check if the numerical result for a squared matrix element at a given phase space point is the same. Example: a quark is replaced by an anti-quark w.r.t. to the other process.
For a set of processes that can be identified by this it is enough to compute one to know them all.
In such a case, the corresponding matrix element squared is calculated only once and then recycled by the other processes.

Colour and spin correlations
The starting point of the Catani-Seymour algorithm is detailed in Eqs. (48) and (49), supplemented with expressions like the one in Eq. (22) for the individual dipole subtraction terms. The latter states that for any given process the Catani-Seymour dipole subtraction term for the real (m + 1)parton correction term consists of the corresponding m-parton matrix element at Born level plus an additional operator that acts on colour and spin space. For the latter, only the limit ǫ → 0 needs to be considered.

• Colour operator:
In all four dipoles, Eq. (25), (31), (35), and (39) colour-correlated tree-amplitudes of the form occur, where i labels the emitter and k the spectator. Denoting the colour indices of the external legs of the tree process explicitly by a j and b j , this can be cast into where T c ab = if acb , if the associated particle is a gluon, and T c ij = t c ij , if the associated particle is a quark. In other words, the colour structure for dipole terms can be generated by adding a gluon connecting the emitter with the spectator as illustrated in Fig. 2. The colour matrix for a dipole term is recomputed after this insertion using the available evaluation tool in AMEGIC++.

• Spin space:
For a quark splitting all spin-matrices are just proportional to δ ss ′ , translating the quark spin to be exactly the same as for the Born-level m-parton matrix element.
For the case of a gluon splitting, however, there are non-trivial correlation matrices. All of them can be cast into the generic form where B andp are functions of the kinematic variables and momenta of the corresponding splitting. Their values are listed in Table 1. dipole type splitting: The structure of the splitting tensor as given in Eq. (63) is very similar to the polarisation sum for massive vector bosons in unitary gauge, except for the factor B and the fact thatp can be timelike or spacelike. This analogy can be used to replace the tensor by a polarisation sum, i.e.
− g µν +p Here the summation index λ runs over four values, +, −, l and s. ξ λ is a sign that cannot be absorbed into the polarisation vectors ǫ λ . For a gauge boson with momentum p µ = p 0 , |˜ p| sin θ cos φ, |˜ p| sin θ sin φ, |˜ p| cos θ , the polarisation vectors are defined as and the sign factors are given by In order to calculate the dipole matrix element, the polarisation vectors of the splitting gluon are then replaced by the ones defined above.

Organisation and process management
To construct all dipole functions necessary to cancel the infrared divergencies of a given parton level real-correction process firstly all pairs of partons have to be determined that might emerge from the splitting of an emitter parton (initial state partons are charge conjugated for this procedure). This might be any quark (or anti-quark) and a gluon, two gluons or a quark and an anti-quark of the same flavour. Secondly, each of those pairs is combined with any possible third parton (acting as spectator) to define all possible dipole functions.
Any individual dipole function is thus specified by: 1. type (the specific combination of initial and final state for emitter and spectator),

the specific flavours involved in the splitting, and
3. the corresponding m-parton matrix element and its emitter and spectator particles.
In order to construct the individual dipole functions, given by the following ingredients are necessary: 1. A rule to map the (m + 1)-parton phase space onto an m-parton phase space.
2. The corresponding splitting function for the dipole. This consists of two parts, a scalar function F (...) of the kinematic variables of the splitting and a spin correlation matrix. As discussed above, for quark splittings the matrix is simply δ ss ′ , for gluon splitting the matrix is represented by an outer product of pseudo-polarisation vectors, which are also functions of the kinematic variables of the splitting.
3. The colour matrix C ′ ij , respecting the extra colour correlation.

Amplitudes
A i of the corresponding m-parton matrix elements. For gluon splitting cases these amplitudes have to be calculated replacing polarisation vectors of the splitting gluon by the pseudo-polarisation vectors introduced above.
The calculation of any dipole function is organized in the class Single DipoleTerm, each instance of this class representing one dipole. This class controls the ingredients for the calculation: Firstly there is a Born-level m-parton matrix element of the original AMEGIC++ implementation, just extended such that it includes the additional colour correlation. Secondly there is a class Dipole Splitting Base that completely organizes the splitting function itself. Specified by the type of the dipole (initial and final states for emitter and spectator) and the type of the splitting (determined by the contributing flavours) it takes care of the mapping between the m + 1-parton and the m-parton phase spaces and of the calculation of the splitting function (including the polarisation vectors to encode the spin correlation).
Above that the class Single Real Correction handles all contributions to an infrared regularized parton level process. This consists firstly of an (m + 1)-parton tree level matrix element in the original AMEGIC++ implementation. Secondly it contains a list of single dipole functions, simply determined by looping over all partons and selecting valid dipole configurations. The classes Single Real Correction and Single Process are derived from a common base class in a way such that the class Process Group can be reused to also organize the infrared regularized parton level process in groups of common features up to all subprocesses contributing to a jet cross section.
Similarly to the case of tree level processes in AMEGIC++, also here a mapping of parton level processes that lead to identical or proportional results can be used to speed up the calculation and save computer resources. To this end, the following automatic identification strategies are implemented: • If two real correction processes can be mapped (using strategies described in section 3.1.2) then also the whole Single Real Correction is mapped.
• For single dipole terms a unique identification algorithm proceeds as follows: Two terms can be mapped if the included m-parton process can be mapped and if the three particle labels (numbering the the external particles of the real correction process) to identify a dipole are identical.
• Many of the born matrix elements within the dipole terms will be identical. However, since different dipoles require different momentum mappings they have to be recalculated. Only the calculation routine can be shared.

Analytical structure of the full result
The starting point of the discussion of the finite pieces of the integrated dipole terms is Eq. (52), where now the phase space integration as well as the summation and integration over the incoming parton flavours and momenta is made explicit. Then, terms inside the m-parton integral come from subtraction terms integrated over the phase space of the extra parton emission and from the collinear counterterm for the general case of a NLO cross section with initial state partons. The terms inside the (m + 1)-parton phase space integral in contrast corresponds to the dipole subtraction bit. Altogether, and including the convolution with parton distribution, the relevant term to be evaluated can thus be cast into The only correlation of the insertion operators I, P, and K with the Born level matrix element is within color space. To be more specific, this implies that only the following structures emerge for all i = j, where i and j may label both final and initial state partons. Since any of the appearing matrix elements with insertion operators can be written as a sum of such structures, the color factors will be skipped in the following and the operators will be treated simply as scalar functions. The terms P and K induce dependences on x, which combined yield result in the structure Here, h(x) and k(x) are regular functions in x and the '+'-distribution is defined by its action on a generic test function a(x) Then the r.h.s. of Eq. (69) can be cast into the form The functions g a,a ′ (x), k a,a ′ (x), and h a,a ′ (1) can be read off the corresponding functions in App. A.
Computationally the most demanding part is the actual Born-level cross section dσ B ab , due to its potentially expensive multi-particle matrix element, which typically suffers from factorial growth with the number of external particles. Thus, the calculation can be significantly accelerated if the expression is rearranged such that dσ B ab has to be computed only once for a single configuration at a given phase space point. This can be achieved by changing the integration variables η to η ′ = xη. After renaming η ′ back to η and reordering the summation over a and a ′ (b and b ′ ) the expression above reads where the G a,b (η) = η 0 dx g a,b (x) are analytically computed. The insertion operator I(ǫ), Eq. (53) is given as a Laurent series in ǫ. For the implementation the interesting part is ∝ ǫ 0 , since the poles must have been analytically extracted before 4 .

Implementation and Organisation
The numerical calculation of the finite contributions from integrated counterterms is organized as Eq. (74) suggests, i.e. the basic unit (class Single Virtual Correction) covers everything that is associated with a specific m-parton cross section. For the actual calculation, basically all colour correlated matrix elements in Eq. (70) are necessary. The contributing amplitudes are, of course, the same for all of them, only the colour matrix is different. Therefore, a generalized version of Single Process is employed that is able to deal with a multitude of colour matrices to calculate all required matrix elements at once. Anything else needed for the calculation of the finite contribution is a long list of rather simple scalar functions and constants. The integration over x is done numerically, i.e. for each set of external momenta x is diced within the corresponding interval.

Phase space integration
Together with the automatic generation of matrix elements AMEGIC++ also generates specific, process-dependent phase-space mappings for efficient integration. There, some a priori knowledge about the integrand is used [60]- [62] together with self-adaptive Monte Carlo integration methods [16,63]- [65]. Here, the general method will briefly be summarized for the case of tree-level processes and its application to the integrals coming with the subtraction method will be discussed.

Importance sampling and Multi-channel integration
The general idea behind importance sampling is to improve the numerical behaviour of an integrand by a change of integration variables, The new variable y is chosen in a way such that f g is a sufficiently smooth function, leading to a reduced error estimate of the integration. Typically, the weight g is chosen as a simplification/approximation of f , such that the integral y = gdx can be analytically solved.
For phase space integrals maps relating vectors of uniformly distributed random numbers {a i } inside the interval [0, 1] to the four-momenta of the external particles of a physical process {p j } are in the center of the sampling process. The weight function in such a case g is then determined by For a single squared amplitude it is easy to determine suitable momentum mappings and weights: invariant masses of the propagators are determined according to the propagator and angles for particle splitting are chosen isotropically; the finial state momenta are then determined out of those variables. For instance, for a massless propagator the invariant s would be generated by with the corresponding weight The constants s max and s min are upper and lower boundaries of the invariant mass, which depend on the overall topology of the phase space point and potential cuts; ν in contrast is an effective exponent for the propagator, subject to choice.
Weight distributions for contributions from several amplitudes can be then combined using the multi-channel method. A total weight function G is defined through where the g k are the weight functions for the single contributions (channels) and the α k are arbitrary coefficients with α k > 0 and k α k = 1. The corresponding momentum mapping is then given by The multi-channel method relies on automatically adapting the coefficients α k such that the variance of the phase space integral is minimized.

Further refinement
The efficiency of the integrator is improved if additionally the self-adaptive VEGAS algorithm [64] is applied on the channels. VEGAS is very efficient in the numerical adaptation to functions, where the peaking behaviour is not too extreme and which are factorisable to a product of one-dimensional functions. This is clearly not given for full matrix elements. However, the structure represented by a single channel fulfills this condition. Thus, in AMEGIC++ VEGAS is used to adapt selected channels to structures that go beyond their rough approximations and which are typically hard to specify analytically or which are a priori unknown.
For each channel VEGAS is used to generate a mapping ξ from uniformly distributed random numbers to a non-uniform distribution, still inside the interval [0, 1], and a corresponding weight v k . To combine this with the multi-channel method the mapping X({a i }) for single channels must meet the requirement to be invertible. The full map reads For a momentum configuration {p j } the weight is therefore given by

Subtracted processes
The subtraction method necessitates the evaluation of two independent integrals, namely integrals over the m-parton and the (m + 1)-parton phase space. In both cases mappings generated for the tree level process of the same dimensionality are used.
For the integration of the (m + 1)-parton phase space soft and collinear regions must be included. In this case the lower limit for the invariant masses of many propagators (e.g. Eq. (78)) must be zero.
To keep the integral over the weight finite the exponent ν must be set to a number smaller than 1. The actual shape of those propagators is hard to specify a priori. It depends on the jet definition and on the balance between the real correction process and the subtraction term (the integrand can be positive or negative). Taken together, however, it seems not unreasonable to assume a small exponent. The VEGAS refinement adapts very good to the actual shape and the final integration efficiency after optimisation has only a weak dependence on the initial value of ν. Since the VEGAS algorithm optimizes on the variance of the integrand it can, to some extend, also deal with the numerical problems related to "missed binning", which will be discussed in the following section.
The m-parton phase space is much simpler. Since most parts of the integrand are proportional to the born matrix element it tends to work very well with this phase space setup.

Cuts and analysis framework for NLO calculations
Triggers and observables for NLO calculations have to be chosen with care. The general strict requirement not to spoil the cancelation of infrared divergencies has already been discussed in section 2.3. Before going into any details concerning cuts, it is important to notice that a rule is mandatory of how cuts act on the different contributions to the NLO cross section. This rule must exist in a m-parton and a (m + 1)-parton version, where the latter needs to satisfy the conditions of infrared safety in degenerate phase space regions. In practical terms, this implies that the (m + 1)-parton version of the cut must allow for exactly one parton to become soft or collinear, while the m-parton version has to omit all singular regions. Second, Eq. (18) requires for the cut of the m + 1 phase space integral to be applied separately to the real correction process (using the m + 1-parton version) and to each dipole term (using the m-parton version, applied on the momenta of the mapped m-parton configuration). In general there might be kinematic configurations, where the real correction process ends up outside the accepted phase space region but some dipole terms do not and vice versa. This leads to the problem of "missed binning": if such a configuration occurs close to a singular region, large contributions result, which do not cancel completely. Ultimately, this leads to large numerical fluctuations, which need to be addressed. This is a common issue for all subtraction methods. So far, the following cuts have been made available in AMEGIC++: • A simple cut for jets is implemented as follows: a suitable jet algorithm (e.g. k T ) [66]- [69] is used to construct jets from the final state partons and their momenta. Then the number of jets above a given p T -cut is counted. A phase space point is valid if this number is greater or equal m.
• Of course also cuts that only act on particles not taking part in strong interactions can be applied. If initial-initial dipoles are present this also has to be done separately for the real correction and for the dipole terms, since the momentum mapping in this case modifies all final state particle momenta. Implemented are cuts on invariant masses, on total or transverse energies, on rapidities or on particle angles w.r.t. the beam.
Sherpa's ANALYSIS-package has been extended to be able to deal with weighted events from the NLO subtraction procedure. For example, and to be more specific, consider the case of a cross section which is differential to some infrared safe quantity F , i.e. a distribution to be binned in a histogram dF . For the m-parton integral no special treatment is mandatory: for a given momentum configuration, dF can directly be evaluated and filled into the corresponding bin. For the real correction and the dipole subtraction functions in the (m + 1)-parton integral, F has to be evaluated for each contribution separately, similar to the phase-space cut. Again, the problem of "missed binnings" appears, if contributions to a single event end up in more than one bin.

Checks of the implementation
In this section a number of tests of the correct implementation of the subtraction algorithm and of the integration routines are described. These tests are mainly technical in nature, results relating to truly physical observables are discussed in the next section, Sec. 5.

Explicit comparisons
Before moving on to technical checks, it is worth stating that a number of direct comparisons of individual terms from the program presented here with those obtained from M. Seymour's Fortran code DISENT have been performed. The latter is a dedicated program to compute NLO cross sections for the deep inelastic scattering processes e − p → e − +jet, e − p → e − +2jets and for electron-positron annihilation to two and three jets. This direct comparison is possible, since DISENT uses exactly the same subtraction formalism, allowing to compare individual terms at given phase space points. All terms listed in the following showed full agreement of the two codes, up to the numerical precision. The comparison included: • Dipole subtraction terms for the real correction: all flavour configurations for dipoles with final state emitters and spectators as well as for dipoles with initial state emitters / final state spectators and final state emitters / initial state spectators have been checked.
• Terms from the finite part of the insertion operator I, cf. Eqs. (49) and (53).
• Terms from the insertion operators K and P for the case of one initial state parton, cf. Eq. (52) and the implemented version Eq. (74).
Furthermore integrated results of the virtual and real parts of the NLO corrections in this subtraction scheme where compared and agreed within statistical errors for all accessible processes.

Test of convergence for the real ME
An obvious first technical check of the overall package consists of testing the convergence behaviour of the dipole subtraction terms close to the singular region. To this end, the m + 1-parton phase space of the regularized real correction part is numerically integrated over. The crucial issue is to ensure that the integrand remains finite over the full phase space, in addition the performance of the integration algorithms deserve some consideration.
Clearly, for the numerical calculation a small phase space region around each singular configuration has to be cut out. Although the dipole terms are expected to become equal to the matrix elements there, technically speaking infinite or very large numbers must be subtracted in this region, leading to large fluctuations and hence to errors due to the limited numerical precision at which the calculation is performed. Therefore a variable α min is introduced, which on the basis of kinematic variables of  corresponding dipole functions, reads as follows: .

(85)
This parameter α serves as a cut-off in such a way that for an externally given parameter α cut kinematic configurations with α min < α cut are omitted.
In Fig. 3 the dependence of the subtracted cross section on α cut for four sets of real correction processes, namely e − e + → 3jets, e − e + → 4jets, e − p → e − + 2jets and pp → 3jets. All types of dipoles and splitting functions contribute to the dipole terms which are necessary to regularize those processes. It is apparent that for α cut ∼ 10 −5 the cross section stabilizes close to its final value.
To study the numerical behaviour near the singularity in more detail, in Fig. 4 the absolute value of the subtracted cross section, binned in intervals of α min is depicted. For all studied processes the contribution to the cross section drops down by at least four orders of magnitude with decreasing α min and confirms the observations for full subtracted cross sections made before.
The strong increase accompanied with statistical errors of 100% or larger for α min values below 10 −9 − 10 −11 signals defects due to the limited numerical precision (double precision ∼ 10 −12 ). One reason is the already mentioned numerical problem when subtracting extreme large and almost equal numbers. Another reason is the precision of the momentum four-vectors itself, because the precision of the external particles residing on their m = 0 mass shell is also limited by the numerical precision. This of course may consequently lead to errors of that order in the matrix element calculation. Thus, Fig. 4 allows to determine best choices for α cut , somewhere between 10 −9 and 10 −11 .

Consistency checks with free parameters
In section 2.6 ways of modifying the subtraction terms without changing the singular behaviour have been discussed. Such modifications can be employed for non-trivial tests of the implementation, since the modifications will affect both, the real part and the virtual part of the NLO cross sections, with their sum remaining constant. In Fig. 5 the total NLO correction for the cross sections of e − e + → 2jets, e − e + → 3jets, e − p → e − + jet and pp → W − → e −ν e and their real and virtual contributions are displayed as functions of the parameter α, as introduced in section 2.6. The fact that the sum remains constant within statistical errors provides a non-trivial confirmation of the correct implementation of the algorithm. It should be noted here that the calculation of the cross section of the processes under consideration invokes all types of dipole functions as well as the most general case of the insertion operators from the integrated dipole terms. By using the same number of phase space points for each integral and comparing statistical errors, it can be seen that this parameter can also be used to optimize the numerical behaviour. Clearly, best results are obtained if the values of the virtual and real contributions are both as small as possible, thus reducing the size of the fluctuations. It should be noted that the error bars in Fig. 5 are given not including the leading order part of the cross section. Since relative errors for the latter can be expected to be much smaller if evaluated for the same number of phase space points, the (relative) statistical error for the full NLO cross section will be significantly reduced.  (57)) for the total cross sections of (a) e − e + → 2jets, (b) e − e + → 3jets, (c) e − p → e − + jet and (d) pp → W − → e −ν e . The results and error bars in the difference plots are determined after calculating 500000 phase space points of the real contribution, which typically dominates the total statistical error.

First physical applications
In this section, some simple applications will demonstrate the performance of the dipole subtraction procedure, as implemented, for the calculation of physically relevant observables. The born matrix elements, dipole subtraction terms to regularize the real correction and corresponding finite terms to be added to the virtual correction were generated automatically by AMEGIC++. The one-loop amplitudes have been explicitly implemented for the considered processes.

Three-jet observables at LEP
To compute three jet cross section at next-to-leading order the one-loop matrix element given in [36] has been implemented. The expression given there is averaged over the direction of incoming momenta, which is sufficient for observables that are not correlated to the beam direction. In Fig. 6 LO and NLO predictions are displayed for observables sensitive to O(α S ). In particular, the event shape observables 1-Thrust, Major, C-parameter and the Durham 3 → 2 jet rate are compared with measurements performed at LEP on the Z 0 -peak by DELPHI [70]. All data are normalized to unity. The normalisation for the calculated cross section, however, is somewhat complicated. This is because in the calculation three-jet events are required in each case, translating into the necessity to apply a phase space cut. On the other hand, the data are more inclusive and also include comparably soft regions, where fixed-order perturbation theory is known to fail and must be supported by resummation techniques. The normalisation for the calculations has thus be chosen such that it agrees with data in the "safe" regions. This exposes the differences between LO and NLO calculations best. As a consequence, the corresponding normalisation factor of both calculations is not identical. From the results of Fig. 6 it can be deduced that for all observables the range described sufficiently well by the calculation is extended for the NLO calculations. For both, the region described by soft physics (left side in all plots) as well as phase space regions populated by additional hard QCD radiation (right side in all event shape plots) the prediction has been improved.

DIS
The one loop matrix element for this process is given by the well known expression where Q 2 = −q 2 > 0 with q the momentum transfer between the electron and the proton. Fig. 7 shows differential cross sections w.r.t. the transverse momentum and rapidity of the scattered electron and the hardest jet at leading and at next-to-leading order. The CM-energy has been taken as √ 10 5 GeV, corresponding to a 50 GeV electron beam and a 500 GeV proton beam. The parton distribution function CTEQ6M [71] has been employed, factorisation and renormalisation scales have been fixed to Q 2 . A phase space cut on the electron of p T > 10GeV has been imposed. The NLO correction for this setup is comparably small, for the total cross section it is of the order of 5% and negative. The ratios of NLO and LO calculation, however, are not constant for all observables. At NLO the cross section rises for increasing momentum transfer, up to a correction of 40% for transverse momenta of electron and jet of the order of 150 GeV.  : Distribution of transverse momentum (left plots) and rapidity, defined in the beam CM frame (right plots), of the scattered electron and of the hardest jet in deep inelastic scattering, calculated at leading order and next-to-leading order. The CM-energy has been taken as √ 10 5 GeV. A phase space cut on the electron (p T > 10GeV ) has been applied. For the rapidity distribution of the first jet a p T > 15GeV has been required. Dashed and dotted lines denote the real and the virtual corrections to the Born cross section, respectively. The lower panels of each plot show the ratio between the leading order and the next-to-leading order results.

W − production at Tevatron
The one-loop virtual contribution to this process can be obtain by crossing relations from Eq. (86) and is given by where now Q 2 =ŝ, the CM energy squared of the incoming partons. Fig. 8 shows cross sections for Tevatron Run II, differential in the rapidity of the W − -boson and the transverse momentum of the electron, respectively. The parton distribution function CTEQ6M [71] has been employed, factorisation and renormalisation scales have been fixed to m 2 W . The total and differential cross sections are in full agreement with predictions obtained using the next-to-leading order parton level generator MCFM [20].

Conclusions and outlook
In this publication a fully automated implementation of the Catani-Seymour dipole formalism in the framework of the matrix element generator AMEGIC++ has been presented. It allows to automatically generate the process-dependent real correction terms for given Born cross sections with massless external particles and the corresponding real subtraction terms. The integration of the subtracted real correction terms is performed automatically with a multi-channel method, giving rise to an appreciable convergence. The implementation has carefully been checked for correctness, invoking consistency checks with free finite terms which may be added to the subtraction terms. Through the explicit inclusion of virtual terms a parton-level calculator is so available.
In the future, the code will be further updated to include massive external particles and to provide a full parton-level generator at NLO.

Acknowledgments
The authors would like to thank M. Seymour for fruitful discussions and numerous comparisons to the code DISENT. Financial support by BMBF and the Marie Curie research training network MCnet (contract number MRTN-CT-2006-035606) is gratefully acknowledged. Furthermore, T. G. would like to thank the Marie Curie Fellowship program for Early Stage Training hosted by CERN for financial support over an extended period of 11 months and the CERN theory group for kind hospitality.
A Insertion operators etc.
In this appendix, the ingredients of the master equation Eq. (52), dσÃ ab (p a , p b ) + dσ C ab (p a , p b , µ 2 F ) = dσ B ab (p a , p b ) × I(ǫ) will be repeated. a and b specify initial state partons, and the sum runs over all accessible a ′ and b ′ occurring in the PDF. The insertion operator I is given by cf. Eq. (53), and again the indices I and J run over all initial and final state partons, while the universal functions V I (ǫ), encoding the singularity structure, merely depend on the flavour of I and read V I (ǫ) = T 2 I 1 ǫ 2 − π 2 3 + γ I 1 ǫ cf. Eq. (54). The individual γ I and K I will be listed in Eqs. (94) and (95).
The factorisation scale dependent terms are proportional to insertion operators P a,a ′ ({p}, xp a , x; µ 2 F ), which read P a,a ′ ({p}, xp a , x; µ 2 The regularized Altarelli-Parisi kernels P ab (x) are listed in Eq. (96).
The factorisation-scheme dependent terms are proportional to the initial-state insertion operators K. For one initial-state hadron only, this operator reads with the functions K aa ′ F.S. (x) andK aa ′ (x) given below, cf. Eqs. (97) and (99), and with the γ i listed in Eq. (94). Note that the subscript "F.S." denotes the factorisation scheme. For two initial state partons, the initial-state insertion operator is given by Eq. (55), with the functionsK aa ′ (x) given in Eq. (98).
The γ I and K I occurring in Eqs. (90) are related to integrals of the Altarelli-Parisi kernels listed below, Eq. (96), and read and respectively. The Altarelli-Parisi kernels emerging in the factorisation-scale dependent terms of Eq. (91) are P qg (x) = Pq g (x) = C F 1 + (1 − x) 2 x P gq (x) = P gq (x) = T R x 2 + (1 − x) 2