text:
# THE TWO CUTOFF PHASE SPACE SLICING METHOD

## I Introduction

Perturbative quantum chromodynamic (QCD) calculations are essential in the effort to describe large momentum transfer hadronic scattering processes. At one time it was sufficient to work at lowest order for the hard scattering subprocesses and utilize the leading-logarithm approximation to treat the higher order gluon radiation and quark-antiquark pair production which give rise to the scale dependence of the parton distribution and fragmentation functions, and to the running of the strong coupling $`\alpha _s`$. As the experimental systematic and statistical errors decreased, the need for increased precision for the theoretical calculations became apparent, leading to the widespread use of next-to-leading-order expressions for the hard scattering subprocesses with the remaining higher order terms being treated in the next-to-leading-logarithm approximation. Early calculations of this type were typically performed with a combination of analytic and numerical integration techniques. The phase space integrations at the parton level were often performed analytically, and the convolutions with the parton distribution or fragmentation functions done numerically. This approach is satisfactory for fully or singly inclusive cross sections, but information is lost about quantities over which the integrations have been performed. Thus, if cuts are to be placed on two or more partons (or hadrons or jets), the calculation must be started anew. Furthermore, for some observables it is difficult to calculate the appropriate Jacobian for the transformation from partonic to hadronic variables. For these reasons it was recognized that Monte Carlo techniques would be useful for such calculations. The Jacobians would be handled by the choice of histogramming variables and several observables could be histogrammed simultaneously. Additionally, it would be simple to define jets and to implement experimental cuts on the four-vectors of the produced partons.

In light of the above observations, a method for performing next-to-leading-logarithm calculations using Monte Carlo techniques was developed bergmann . Two cutoff parameters serve to separate the regions of phase space containing the soft and collinear singularities from the non-singular regions; nowadays this is referred to as the phase-space slicing technique.

The usefulness and generality of the method may be appreciated by considering the many physics processes to which it has been applied. The basic core of the method was first developed to study QCD corrections to dihadron production bergmann . It has subsequently been applied to direct jet photoproduction boo1 , hadronic photon–jet boo2 , direct photon boo3 , W breno , ZZ oo , WW ohnemus1 , WZ ohnemus2 , two photon boo\_2photon ; cg1 , Z$`\gamma `$ ohnemus3 , and W–Higgs os\_whiggs ; bbo\_whiggs  production, nonstandard three vector boson couplings in W$`\gamma `$ bho1 , WZ bho2  and WW bho3  production, hadronic photon–heavy quark production bbg1 ; bbg2 , jet photoproduction ho , quantum electrodynamic (QED) corrections to hadronic Z production bks , QCD corrections to slepton pair production bhr , electroweak corrections to W production bkw , single-top-quark production harris , and dihadron production Owens:2001rr .

Despite this usefulness, a full and detailed description of the method does not exist in the literature. Here we provide such a description. Naturally, as the method was applied to the above physics processes, refinements were made. We therefore take this opportunity to modernize and systematize the presentation relative to that given in bergmann , and show details, which have not previously been published. Searches for signals of new physics often rely on next-to-leading-order Monte Carlo-based calculations. It is anticipated that the details provided here will prove to be helpful for anyone wanting to apply the method to additional processes.

In the course of a next-to-leading-order calculation ultra-violet singularities show up in loop integrals where the momenta go to infinity. They are removed through the process of renormalization. (See, for example, Refs. smith ; sterman  for a discussion.) Soft (infrared) divergences arise if the theory includes a massless field like the photon in QED or the gluon in QCD. They are encountered in both loop and phase space integrals and are found in the low energy region where the integration momenta go to zero. The soft singularities cancel between the virtual and bremsstrahlung processes kln . If the massless field couples to another massless field, or to itself, a collinear (mass) singularity may occur in both loop and phase space integrals. Final state mass singularities cancel when summed over degenerate (experimentally indistinguishable) final states according to the theorem of Kinoshita-Lee-Nauenberg kln . For tagged hadrons there is no final state sum, and the associated mass singularities are factorized into fragmentation functions. Similarly, initial state singularities do not cancel because there is typically no sum over degenerate states; they are removed by factorization css ; bodwin .

The goal of the practitioner of next-to-leading-order calculations is to organize the soft and collinear singularity cancellations described above without loss of information in terms of observable quantities. The phase space slicing method provides a relatively simple and robust method to do this. Several other methods for handling the organization of the cancellations exist in the literature and have been used to study a wide variety of high energy processes, including some of those listed above. The phase space slicing method of one cutoff first developed in ycut1 ; ycut2  divides the phase space according to $`s_{ij}=(p_i+p_j)^2>ys_{12}`$ where $`p_i`$ and $`p_j`$ label the momenta of partons $`i`$ and $`j`$, and $`y`$ is a small dimensionless parameter. Another variant for jets gg ; ggk  and hadrons and heavy quarks kl  partitions phase space according to $`s_{ij}>s_{\mathrm{min}}`$ where $`s_{\mathrm{min}}`$ is a small dimension-full parameter. It is also possible to engineer the singularity cancellation using plus distributions, commonly referred to as the subtraction method, which has been applied to jets ert ; kn ; ks ; eks1 ; eks2 ; Frixione:1995ms ; Frixione:1997np ; Nagy:1996bz  and heavy quark final states mnr ; fmnr ; hs ; ho2 . The subtraction method taken together with factorization formulae that interpolate between the soft and collinear approximations to the matrix elements is known as the dipole method dipole . The dipole method was originally developed for jet and light hadron cross sections where it has seen extensive application Nagy:1997yn ; Weinzierl:1999yf ; Nagy:2001xb ; Campbell:1999ah ; Ellis:1998fv ; Beenakker:1999xh ; Maina:1996ep ; Catani:1999nf ; Kramer:1999bf . It has recently been extended to handle massive fermions and partons ditt ; Phaf:2001gc ; Catani:2002hc  and is finding many applications harris ; Denner:2000bj ; Dittmaier:2001ay ; Dittmaier:2000tc  in that domain as well. A brief comparison of the slicing and subtraction methods is given in Appendix A. Properly implemented, all methods should give identical physics predictions.

This paper proceeds as follows. In Section II we give details of the phase space slicing method with two cutoffs. We examine the soft and collinear regions of phase space and see how to arrive at a finite cross section. In Section III the process of electron-positron annihilation into quarks is studied for massive, massless, and tagged final states. After that, the examples of lepton pair production and single particle production in hadronic collisions are given. We conclude in Section IV. As mentioned above, Appendix A contains a brief comparison with the subtraction method. Appendix B contains angular integrals useful in the soft analysis. A discussion of terms that vanish like the ratio of the two cutoffs is given in Appendix C. Appendix D explains how to improve numerical convergence.

## II The Method

This section contains the main derivations for jet, fragmentation, and heavy quark final states, as well as a discussion of initial state mass factorization. Before getting into too many details, it will be helpful to outline the procedure first. The typical calculation involves lowest order two-to-two subprocesses which have two-body final states and higher order two-to-three subprocesses which lead to both two- and three-body final states. In addition, the one-loop virtual corrections also contribute to the two-body final states.

We begin by decomposing the three-body phase space used to calculate the two-to-three contribution to the partonic cross section into two regions which we call soft, S, and hard, H, by writing 111Implicit in the cross section is a measurement function which serves to implement the jet algorithm and/or define the experimentally visible portion of phase space (the cuts). For the cancellation of singularities to take place, the measurement function is required to be infrared-safesw .

$$\sigma =\frac{1}{2\mathrm{\Phi }}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3=\frac{1}{2\mathrm{\Phi }}_\mathrm{S}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3+\frac{1}{2\mathrm{\Phi }}_\mathrm{H}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3,$$
(1)
where $`\mathrm{\Phi }=\lambda ^{1/2}(s,m_1^2,m_2^2)`$ is the usual flux factor which depends on the partonic center-of-momentum energy squared $`s`$ and the incident particle masses $`m_1`$ and $`m_2`$, $`\overline{}|M_3|^2`$ is the two-to-three body squared matrix element averaged (summed) over initial (final) degrees of freedom, and $`d\mathrm{\Gamma }_3`$ is the three-body phase space. The partitioning of phase space into S and H depends on a parameter $`\delta _s`$ in a manner to be described below. Within S the double pole (eikonal) approximation to the matrix elements is made and then analytically integrated over the unobserved degrees of freedom in $`n`$ space-time dimensions. The result, depending on the masses of the partons involved, may contain double and/or single poles in $`n4`$, and accompanying double and/or single logarithms in the soft cutoff $`\delta _s`$. We always work in the approximation where the cutoffs are small, so terms of order $`\delta _s`$ may be neglected. Just how small is needed will be studied below.

Next, if there are collinear singularities present, the hard region is further decomposed into collinear C, and non-collinear $`\overline{\mathrm{C}}`$, regions as follows:

$$\frac{1}{2\mathrm{\Phi }}_\mathrm{H}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3=\frac{1}{2\mathrm{\Phi }}_{\mathrm{HC}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3+\frac{1}{2\mathrm{\Phi }}_{\mathrm{H}\overline{\mathrm{C}}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3.$$
(2)
This partitioning depends on a second cutoff $`\delta _c`$. Within HC the leading collinear pole approximation to the squared matrix element is made. As explained below, exact collinear kinematics may be used to define the integration domain of HC when $`\delta _c\delta _s`$. The integrations over the unobserved degrees of freedom are performed analytically in $`n`$ space-time dimensions giving a factorized result where single poles in $`n4`$, and single logarithms in both cutoffs $`\delta _c`$ and $`\delta _s`$, multiply splitting functions and lower-order squared matrix elements. This is in the approximation where terms of order $`\delta _c`$ and $`\delta _s`$ are neglected.

The cancellation of poles in $`n4`$ is based on the Kinoshita-Lee-Nauenberg theorem kln , or mass factorization css ; bodwin , depending on the situation. For experimentally degenerate final states, the soft and final state hard collinear singularities cancel upon addition of the interference of the leading order diagrams with the renormalized one-loop virtual diagrams. The remaining initial state collinear singularities are factorized and absorbed into the parton distribution functions. The result is finite in $`n=4`$ dimensions, but depends logarithmically on the cutoffs. For tagged hadrons there is no final state sum, and the associated mass singularities are factorized into fragmentation functions.

The integration over the hard non-collinear $`\mathrm{H}\overline{\mathrm{C}}`$ portion of the phase space is performed using standard Monte Carlo techniques. The result is finite by construction and the expressions may be evaluated in four dimensions. At this stage, the calculation yields a set of two-body weights which have explicit logarithmic dependence on the two cutoffs and the three-body weights for which a logarithmic dependence on the cutoffs develops as the Monte Carlo integration is performed. When all of the contributions are combined at the histogramming stage, the cutoff dependence cancels for suitably defined infrared-safe observables. In the following subsections we look at each of these steps in detail.

### II.1 Soft

In this subsection we describe the procedure for extracting soft gluon singularities. When one of the gluons is soft, the phase space is greatly simplified and the eikonal (double pole) approximation of the matrix element is valid gy ; cs ; bcm ; fmr ; cg . The cross section is simple enough to be analytically integrated over the unobserved degrees of freedom in $`n`$ space-time dimensions. The required integrals are well known hs ; been ; willy . The result, depending on the masses of the partons involved, contains double and/or single poles in $`n4`$, and accompanying double and/or single logarithms in the soft cutoff $`\delta _s`$. We work in the approximation where terms of order $`\delta _s`$ are neglected.

Generically, we begin by writing the two-to-three body contribution to the partonic cross section

$$\sigma =\frac{1}{2\mathrm{\Phi }}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3,$$
(3)
as the sum of soft, $`\mathrm{S}`$, and hard, $`\mathrm{H}`$, terms

$$\sigma =\sigma _\mathrm{S}+\sigma _\mathrm{H},$$
(4)
where

$$\sigma _\mathrm{S}=\frac{1}{2\mathrm{\Phi }}_\mathrm{S}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3,$$
(5)
and

$$\sigma _\mathrm{H}=\frac{1}{2\mathrm{\Phi }}_\mathrm{H}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3.$$
(6)
In this section we examine $`\sigma _\mathrm{S}`$ in detail. Further evaluation of $`\sigma _\mathrm{H}`$ is deferred to the following subsections.

Let the particles in the scattering be labeled by their four-momenta $`p_1+p_2=p_3+p_4+p_5`$ and define the Mandelstam invariants $`s_{ij}=(p_i+p_j)^2`$ and $`t_{ij}=(p_ip_j)^2`$. Consider the case when parton 5 is a soft gluon. The soft region $`\mathrm{S}`$ is defined in terms of the gluon energy $`E_5`$ in the $`p_1+p_2`$ rest frame by $`0E_5\delta _s\sqrt{s_{12}}/2`$. The hard region $`\mathrm{H}`$ is the complement: $`E_5>\delta _s\sqrt{s_{12}}/2`$. The gluon energy can be calculated from the other invariants in the problem as follows. Start with $`p_1+p_2p_5=p_3+p_4`$ which, after squaring both sides, yields $`(p_1+p_2)^22p_5(p_1+p_2)=(p_3+p_4)^2`$. In the $`p_1+p_2`$ rest frame $`p_1+p_2=\sqrt{s_{12}}(1,0,0,0)`$, so $`s_{12}2E_5\sqrt{s_{12}}=s_{34}`$. Solving for the gluon energy gives

$$E_5=\frac{s_{12}s_{34}}{2\sqrt{s_{12}}}.$$
(7)
This expression for $`E_5`$ and the definition of the soft region are independent of the masses of the other particles in the reaction. Now that we have defined the boundaries of the soft portion of phase space, we examine the approximations that can be made.

The three-body phase space in $`n`$ dimensions is given by

$$d\mathrm{\Gamma }_3=\frac{d^{n1}p_3}{2p_3^0(2\pi )^{n1}}\frac{d^{n1}p_4}{2p_4^0(2\pi )^{n1}}\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}(2\pi )^n\delta ^n(p_1+p_2p_3p_4p_5).$$
(8)
The divergence in the integral of the matrix element over phase space will be at worst logarithmic. Therefore, up to corrections of $`𝒪(\delta _s)`$, we can set $`p_5^\mu =0`$ in the delta function and regroup terms, yielding

$$d\mathrm{\Gamma }_3|_{\mathrm{soft}}=\left[\frac{d^{n1}p_3}{2p_3^0(2\pi )^{n1}}\frac{d^{n1}p_4}{2p_4^0(2\pi )^{n1}}(2\pi )^n\delta ^n(p_1+p_2p_3p_4)\right]\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}.$$
(9)
This is simply

$$d\mathrm{\Gamma }_3|_{\mathrm{soft}}=d\mathrm{\Gamma }_2\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}},$$
(10)
where

$$d\mathrm{\Gamma }_2=\frac{d^{n1}p_3}{2p_3^0(2\pi )^{n1}}\frac{d^{n1}p_4}{2p_4^0(2\pi )^{n1}}(2\pi )^n\delta ^n(p_1+p_2p_3p_4),$$
(11)
is the two-body phase space of partons 3 and 4. Likewise, up to corrections of $`𝒪(\delta _s)`$, we can parameterize the gluon’s $`n`$momentum in the $`p_1+p_2`$ rest frame as

$$p_5=E_5(1,\mathrm{},\mathrm{sin}\theta _1\mathrm{sin}\theta _2,\mathrm{sin}\theta _1\mathrm{cos}\theta _2,\mathrm{cos}\theta _1),$$
(12)
where the dots indicate the $`n4`$ unspecified momentum components which may be trivially integrated over using

$`d^{n1}p_5`$ $`=`$ $`d|\stackrel{}{p_5}||\stackrel{}{p_5}|^{n2}d\mathrm{\Omega }_{n2}`$ (13)
$`=`$ $`dE_5E_5^{n2}\mathrm{sin}^{n3}\theta _1d\theta _1\mathrm{sin}^{n4}\theta _2d\theta _2\mathrm{\Omega }_{n4}.`$
The angular volume element

$$\mathrm{\Omega }_{n4}=\frac{2\pi ^{(n3)/2}}{\mathrm{\Gamma }[(n3)/2]},$$
(14)
may be rewritten using

$$\mathrm{\Gamma }[(n3)/2]=\sqrt{\pi }2^{2ϵ}\frac{\mathrm{\Gamma }(12ϵ)}{\mathrm{\Gamma }(1ϵ)},$$
(15)
where we have set $`n=42ϵ`$. The final result for the phase space volume approximated in the soft region is

$$d\mathrm{\Gamma }_3|_{\mathrm{soft}}=d\mathrm{\Gamma }_2\left[\left(\frac{4\pi }{s_{12}}\right)^ϵ\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\frac{1}{2(2\pi )^2}\right]dS,$$
(16)
with

$$dS=\frac{1}{\pi }\left(\frac{4}{s_{12}}\right)^ϵ_0^{\delta _s\sqrt{s_{12}}/2}𝑑E_5E_5^{12ϵ}\mathrm{sin}^{12ϵ}\theta _1d\theta _1\mathrm{sin}^{2ϵ}\theta _2d\theta _2.$$
(17)
Once we have the corresponding soft approximation to the matrix element this integral can be performed, yielding the advertised singularity structure. But before proceeding, we pause to note that choosing 5 to be the soft gluon is not special. A similar analysis holds when parton 3 or 4 is a gluon. At this order in perturbation theory, only one gluon may be soft at a time.

Soft photon emission in QED is characterized by the factorization of an eikonal current from the scattering amplitude. The structure of multiple soft photon emission has been studied by Grammer and Yennie gy . In QCD the process is different because gluons carry color charge and can therefore radiate during the scattering. Fortunately, when QCD amplitudes are decomposed into color sub-amplitudes they enjoy the same factorization properties as QED amplitudes cs ; bcm ; fmr ; cg .

Let parton 5 be the soft gluon and take it to have color index $`a`$ ($`=1,\mathrm{},N^21`$) and Lorentz index $`\mu `$. The matrix element factorizes as

$$M_3^a|_{\mathrm{soft}}g\mu _r^ϵ\epsilon ^\mu (p_5)𝐉_\mu ^a(p_5)𝐌_2.$$
(18)
The mass dimensions of the strong coupling have been isolated into the parameter $`\mu _r`$ which is identified with the renormalization scale, leaving the dimensionless coupling $`g`$. The soft gluon’s polarization vector, denoted by $`\epsilon ^\mu (p_5)`$, is Lorentz contracted with the non-abelian eikonal current given by

$$𝐉_\mu ^a(p_5)=\underset{f=1}{\overset{4}{}}𝐓_f^a\frac{p_f^\mu }{p_fp_5},$$
(19)
which itself is color contracted with the color sub-amplitude $`𝐌_2`$. The sum corresponds to the soft gluon being emitted from each external line in turn. The $`SU(N)`$ color charge associated with the emitting parton $`f`$ is denoted by $`𝐓_f`$. Squaring and summing Eq. (18) over the soft gluon polarizations gives

$$|M_3|^2|_{\mathrm{soft}}g^2\mu _r^{2ϵ}\underset{f,f^{}=1}{\overset{4}{}}\frac{p_fp_f^{}}{p_fp_5p_f^{}p_5}M_{ff^{}}^0,$$
(20)
where

$$M_{ff^{}}^0=(𝐓_f^a𝐌_2)(𝐓_f^{}^a𝐌_2)^{^{}}=\left[M_{c_1\mathrm{}b_f\mathrm{}b_f^{}\mathrm{}c_4}\right]^{}T_{b_fd_f}^aT_{b_f^{}d_f^{}}^aM_{c_1\mathrm{}d_f\mathrm{}d_f^{}\mathrm{}c_4},$$
(21)
is the square of the color connected Born amplitude. If the emitting parton is a final state quark or initial state anti-quark, the color charge is in the fundamental representation $`(T^a)_{ij}=t_{ij}^a`$ ($`i,j=1,\mathrm{},N`$). For a final state anti-quark or initial state quark $`(T^a)_{ij}=t_{ij}^a`$. If the emitting parton is a gluon, the color charge is in the adjoint representation $`(T^a)_{bc}=if_{abc}`$.

Substituting (16) and (20) into (5) gives

$$d\sigma _\mathrm{S}=\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\underset{f,f^{}=1}{\overset{4}{}}d\sigma _{ff^{}}^0\frac{p_fp_f^{}}{p_fp_5p_f^{}p_5}𝑑S,$$
(22)
where

$$d\sigma _{ff^{}}^0=\frac{1}{2\mathrm{\Phi }}\overline{}M_{ff^{}}^0d\mathrm{\Gamma }_2.$$
(23)
The integration over the eikonal factors depends on the masses of the particles in the reaction. We leave the integrals to be performed on a case-by-case basis, although all possible mass combinations may be worked out. Specific applications of Eq. (22) are given below.

From a practical point of view, the presentation just made requires only the calculation of the eikonal factors and the color connected Born amplitudes in $`n=42ϵ`$ dimensions. There is no need to evaluate the full two-to-three body matrix element in $`n`$ dimensions. Depending on the complexity of the process, this may be a simplification. However, if the full two-to-three body matrix element in $`n=42ϵ`$ dimensions is at hand, setting $`p_5^\mu =0`$ everywhere in the numerator and retaining only the leading singular terms as $`p_5^\mu 0`$ will reproduce Eq. (20) directly, this being known as the double pole approximation bergmann .

### II.2 Collinear

We now return to the further evaluation of the hard portion of the cross section which was separated out in Sec. II.1. The phase space is greatly simplified in the limit where two of the partons are collinear. In the same limit, the leading pole approximation of the matrix element is valid. The cross section is simple enough to be analytically integrated over the unobserved degrees of freedom in $`n`$ space-time dimensions. The result contains single poles in $`n4`$, and accompanying logarithms of the soft $`\delta _s`$ and collinear $`\delta _c`$ cutoffs. Terms of order $`\delta _s`$ and $`\delta _c`$ are neglected.

To this end we further decompose $`\sigma _\mathrm{H}`$ given in Eq. (6) into a sum of hard-collinear $`\mathrm{HC}`$ and hard-non-collinear $`\mathrm{H}\overline{\mathrm{C}}`$ terms

$$\sigma _\mathrm{H}=\sigma _{\mathrm{HC}}+\sigma _{\mathrm{H}\overline{\mathrm{C}}},$$
(24)
with

$$\sigma _{\mathrm{HC}}=\frac{1}{2\mathrm{\Phi }}_{\mathrm{HC}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3,$$
(25)
and

$$\sigma _{\mathrm{H}\overline{\mathrm{C}}}=\frac{1}{2\mathrm{\Phi }}_{\mathrm{H}\overline{\mathrm{C}}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3.$$
(26)
The $`\mathrm{HC}`$ regions of phase space are those where any invariant ($`s_{ij}`$ or $`t_{ij}`$) becomes smaller in magnitude than $`\delta _cs_{12}`$, the collinear condition, while at the same time all gluons remain hard. The complementary $`\mathrm{H}\overline{\mathrm{C}}`$ pieces are finite and may be evaluated numerically in four dimensions using standard Monte-Carlo techniques vegas .

The piece containing the collinear singularities, $`\sigma _{\mathrm{HC}}`$, is treated according to whether the singularities are initial or final state in origin. For the former, factorization provides the formalism for removing the singularities. In the latter case, we distinguish between experimentally degenerate and tagged final states, and rely on either the Kinoshita-Lee-Nauenberg theorem or factorization to dispose of the singularities. We discuss the final state cases first, then return to the initial state.

#### II.2.1 Indistinguishable Final States

Consider the case when there is a sum over experimentally degenerate final states, such as a jet or total cross section. Let partons 4 and 5 be massless and collinear to each other, $`0s_{45}\delta _cs_{12}`$. If we define $`p_{45}=p_4+p_5`$, then for fixed $`p_5`$ we have $`d^{n1}p_{45}=d^{n1}p_4`$. The three body phase space Eq. (8) may be written as

$$d\mathrm{\Gamma }_3|_{\mathrm{coll}}=\left[\frac{d^{n1}p_3}{2p_3^0(2\pi )^{n1}}\frac{d^{n1}p_{45}}{2p_{45}^0(2\pi )^{n1}}(2\pi )^n\delta ^n(p_1+p_2p_3p_{45})\right]\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}\frac{p_{45}^0}{p_4^0}.$$
(27)
This is simply

$$d\mathrm{\Gamma }_3|_{\mathrm{coll}}=d\mathrm{\Gamma }_2\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}\frac{p_{45}^0}{p_4^0},$$
(28)
where $`d\mathrm{\Gamma }_2`$ is the two-body phase space of the particles 3 and 45. In the collinear limit ($`p_t0`$ with $`z`$ fixed) we can write

$`p_{45}`$ $`=`$ $`(P,0,0,P)`$
$`p_4`$ $``$ $`(zP+{\displaystyle \frac{p_t^2}{2zP}},\stackrel{}{p}_t,zP)`$
$`p_5`$ $``$ $`((1z)P+{\displaystyle \frac{p_t^2}{2(1z)P}},\stackrel{}{p}_t,(1z)P).`$ (29)
Then $`p_4+p_5=p_{45}+𝒪(p_t^2)`$ and

$$s_{45}=2p_4p_5\frac{p_t^2}{z(1z)}.$$
(30)
Now using $`d^{n2}p=dpp^{n3}d\mathrm{\Omega }_{n3}`$ and $`z(1z)ds_{45}=dp_t^2`$ we find

$$\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}\frac{p_{45}^0}{p_4^0}=\frac{(4\pi )^ϵ}{16\pi ^2\mathrm{\Gamma }(1ϵ)}dzds_{45}[s_{45}z(1z)]^ϵ.$$
(31)
The corresponding approximation to the matrix element is obtained by imposing collinear kinematics on the portion of the two-to-three matrix element proportional to the leading collinear singularity. This is known as the leading pole or collinear approximation. As a consequence of the factorization theorems css ; bodwin , the squared matrix element factors into the product of a splitting kernel and a leading order squared matrix element bergmann ; ert ; ks ; mnr ; ap .

As above, let partons 4 and 5 be collinear. The matrix element factorizes as

$$\overline{}|M_3(1+23+4+5)|^2\overline{}|M_2(1+23+4^{})|^2P_{44^{}}(z,ϵ)g^2\mu _r^{2ϵ}\frac{2}{s_{45}},$$
(32)
where the $`P_{ij}(z,ϵ)`$ are the unregulated ($`z<1`$) splitting functions calculated in $`n=42ϵ`$ dimensions related to the usual Altarelli-Parisi splitting kernels ap . We label as $`4^{}`$ the parton which splits into the 45 collinear pair. Generally, Eq. (32) contains an additional term that vanishes after integration over the azimuthal angles in $`n`$ dimensions ks ; mnr . Such a term does not contribute to our result.

Substituting Eqs. (28), (31), and (32) into Eq. (25) gives

$`d\sigma _{\mathrm{HC}}^{1+23+4+5}`$ $`=`$ $`d\sigma _0^{1+23+4^{}}\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]`$ (33)
$`\times `$ $`{\displaystyle _0^{\delta _cs_{12}}}{\displaystyle \frac{ds_{45}}{s_{45}}}\left({\displaystyle \frac{s_{45}}{s_{12}}}\right)^ϵ{\displaystyle 𝑑z[z(1z)]^ϵP_{44^{}}(z,ϵ)},`$
where we have used

$$\frac{1}{\mathrm{\Gamma }(1ϵ)}=\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}+𝒪(ϵ^2).$$
(34)
The collinear condition ($`0s_{45}\delta _cs_{12}`$) sets the $`s_{45}`$ integration limits. The hard condition sets the $`z`$ integration limits which depend additionally on the splitting function involved, and also on the mass of the parton recoiling against the 45 pair (3 in this case). This later dependence enters through the threshold condition, as discussed below.

First take parton 3 to be massless. For $`qqg`$ splitting the hard region is defined by $`\delta _s\sqrt{s_{12}}/2E_5\sqrt{s_{12}}/2`$, assuming 4 labels $`q`$ and 5 labels $`g`$. From Eq. (29) we have $`s_{34}=(p_3+p_4)^2=2p_3p_4(2p_3p_{45})z`$ and $`s_{12}=(p_3+p_{45})^2=s_{45}+2p_3p_{45}2p_3p_{45}`$ which together yield $`s_{34}zs_{12}`$. Using Eq. (7) the hard condition becomes

$$0z1\delta _s.$$
(35)
For the $`ggg`$ splitting it is required that both gluons be hard, i.e., $`E_4`$ and $`E_5\delta _s\sqrt{s_{12}}/2`$. $`z`$ then satisfies the relation $`\delta _sz1\delta _s`$. For the $`gq\overline{q}`$ splitting there are no soft singularities, so $`0z1`$ may be taken. In all of these cases the $`z`$ integration limits are independent of $`s_{45}`$ by virtue of the approximation $`s_{45}=0`$ implicit in $`s_{12}2p_3p_{45}`$ (this point is discussed further at the end of Sec. III.2 and in Appendix C). The outermost integration over $`s_{45}`$ may therefore be performed giving the result

$$d\sigma _{\mathrm{HC}}^{1+23+4+5}=d\sigma _0^{1+23+4^{}}\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{1}{ϵ}\right)\delta _c^ϵ𝑑zz^ϵ(1z)^ϵP_{44^{}}(z,ϵ).$$
(36)
For $`z<1`$ the $`n`$dimensional unregulated splitting functions may be written as $`P_{ij}(z,ϵ)=P_{ij}(z)+ϵP_{ij}^{}(z)`$. Explicitly,

$`P_{qq}(z)`$ $`=`$ $`C_F{\displaystyle \frac{1+z^2}{1z}}`$ (37)
$`P_{qq}^{}(z)`$ $`=`$ $`C_F(1z)`$ (38)
$`P_{gq}(z)`$ $`=`$ $`C_F{\displaystyle \frac{1+(1z)^2}{z}}`$ (39)
$`P_{gq}^{}(z)`$ $`=`$ $`C_Fz`$ (40)
$`P_{gg}(z)`$ $`=`$ $`2N\left[{\displaystyle \frac{z}{1z}}+{\displaystyle \frac{1z}{z}}+z(1z)\right]`$ (41)
$`P_{gg}^{}(z)`$ $`=`$ $`0`$ (42)
$`P_{qg}(z)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[z^2+(1z)^2\right]`$ (43)
$`P_{qg}^{}(z)`$ $`=`$ $`z(1z),`$ (44)
where $`N=3`$ and $`C_F=(N^21)/2N=4/3`$ for QCD. Expanding the integrand in Eq. (36) to $`𝒪(ϵ)`$ and integrating over $`z`$ yields the final state hard-collinear terms

$$d\sigma _{HC}^{1+23+4+5}=d\sigma _0^{1+23+4^{}}\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_1^{4^{}45}}{ϵ}+A_0^{4^{}45}\right),$$
(45)
where

$`A_1^{qqg}`$ $`=`$ $`C_F\left(3/2+2\mathrm{ln}\delta _s\right)`$ (46)
$`A_0^{qqg}`$ $`=`$ $`C_F\left[7/2\pi ^2/3\mathrm{ln}^2\delta _s\mathrm{ln}\delta _c\left(3/2+2\mathrm{ln}\delta _s\right)\right]`$ (47)
$`A_1^{gq\overline{q}}`$ $`=`$ $`n_f/3`$ (48)
$`A_0^{gq\overline{q}}`$ $`=`$ $`n_f/3\left(\mathrm{ln}\delta _c5/3\right)`$ (49)
$`A_1^{ggg}`$ $`=`$ $`N\left(11/6+2\mathrm{ln}\delta _s\right)`$ (50)
$`A_0^{ggg}`$ $`=`$ $`N\left[67/18\pi ^2/3\mathrm{ln}^2\delta _s\mathrm{ln}\delta _c\left(11/6+2\mathrm{ln}\delta _s\right)\right].`$ (51)
where $`n_f`$ denotes the number of active flavors.

When the mass of the parton recoiling against the 45 pair, $`m_3`$, is retained, the factorization of the phase space and matrix element is unaffected. Likewise, the collinear condition $`0s_{45}\delta _cs_{12}`$ remains unchanged. It is only the boundries of the hard region that are modified. The full kinematic range of the invariant $`s_{34}`$ is $`m_3^2s_{34}s_{12}`$. The threshold for producing a particle of mass $`m_3`$ sets the lower limit. In terms of $`s_{34}`$ the hard region is $`m_3^2s_{34}(1\delta _s)s_{12}`$. This implies the hard condition Eq. (35) becomes $`0z1\delta _s/(1m_3^2/s_{12})`$. A similar analysis follows for the $`ggg`$ splitting case wherein $`\delta _s/(1m_3^2/s_{12})z1\delta _s/(1m_3^2/s_{12})`$. We therefore immediately see that Eqs. (46)–(51) are valid with the replacement $`\delta _s\delta _s/(1m_3^2/s_{12})`$.

#### II.2.2 Tagged Final States

Next, consider a process where a particular type of hadron is identified in the final state. This necessitates the introduction of a fragmentation function $`D_{h/c}(z)`$ which gives the probability density for finding a hadron $`h`$ which carries a fraction $`z`$ of the momentum of the parent parton $`c`$. Consider the case where parton 4 fragments into a hadron $`h`$, for which the lowest order cross section is

$$d\sigma _0^{1+23+h}=d\sigma _0^{1+23+4}D_{h/4}(z)dz.$$
(52)
The hard collinear cross section expression in Eq. (36) becomes

$`d\sigma _{\mathrm{HC}}^{1+23+h+5}`$ $`=`$ $`d\sigma _0^{1+23+4^{}}\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]\left({\displaystyle \frac{1}{ϵ}}\right)\delta _c^ϵ`$ (53)
$`\times `$ $`{\displaystyle 𝑑yy^ϵ(1y)^ϵD_{h/4}(x)𝑑xP_{44^{}}(y,ϵ)\delta (xyz)𝑑z}.`$
The delta function insures that the hadron $`h`$ carries a momentum fraction $`z`$ of the parent parton’s momentum (parton $`4^{}`$ in this example). Here there is a splitting $`4^{}45`$ followed by parton 4 fragmenting to hadron $`h`$. When all possible $`23`$ subprocesses are considered, there will be several contributions of this same form, corresponding to a sum over the parton 4. For example, if $`4^{}`$ is a gluon, there can be $`ggg`$ followed by $`gh`$ or $`gq\overline{q}`$ followed by $`qh`$ or $`\overline{q}h`$. Similarly, if $`4^{}`$ is a quark $`q`$, there can be $`qqg`$ followed by $`qh`$ or by $`gh`$. Furthermore, the limits of integration on $`y`$ depend on the splitting function as in the case discussed in the previous subsection.

The collinear singularity, evidenced by the pole in $`ϵ`$ in Eq. (53), must be factorized and absorbed into the bare fragmentation function. To do this, we introduce a scale dependent parton fragmentation function

$$D_{h/c}(z,M_f)=D_{h/c}(z)+\left(\frac{1}{ϵ}\right)\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{M_f^2}\right)^ϵ\right]_z^1\frac{dy}{y}D_{h/c^{}}(z/y)P_{c^{}c}^+(y).$$
(54)
In this expression there is an implied sum over the index $`c^{}`$ corresponding to the sum over the different fragmentation possibilities referred to above. The final state factorization scale has been denoted by $`M_f`$. Notice, too, that the integration over $`y`$ extends from $`z`$ to 1. This form for the scale dependent fragmentation function corresponds to the $`\overline{\mathrm{MS}}`$ convention. The regulated ($`x1`$) splitting functions ap  are given by

$`P_{qq}^+(x)`$ $`=`$ $`C_F\left[{\displaystyle \frac{1+x^2}{(1x)_+}}+{\displaystyle \frac{3}{2}}\delta (1x)\right]`$ (55)
$`P_{gq}^+(x)`$ $`=`$ $`C_F\left[{\displaystyle \frac{1+(1x)^2}{x}}\right]`$ (56)
$`P_{gg}^+(x)`$ $`=`$ $`2N\left[{\displaystyle \frac{x}{(1x)_+}}+{\displaystyle \frac{1x}{x}}+x(1x)\right]+\left({\displaystyle \frac{11}{6}}N{\displaystyle \frac{1}{3}}n_f\right)\delta (1x)`$ (57)
$`P_{qg}^+(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[x^2+(1x)^2\right].`$ (58)
Next, we rewrite the bare fragmentation function in Eq. (52) in terms of the scale dependent expression given above, yielding to $`𝒪(\alpha _s)`$

$`d\sigma _0^{1+23+h}`$ $`=`$ $`d\sigma _0^{1+23+4}`$
$`\times `$ $`\left\{D_{h/4}(z,M_f)+\left({\displaystyle \frac{1}{ϵ}}\right)\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{M_f^2}}\right)^ϵ\right]{\displaystyle _z^1}{\displaystyle \frac{dy}{y}}D_{h/c^{}}(z/y)P_{c^{}4}^+(y)\right\}dz.`$
The second term is sometimes refered to as the mass factorization counterterm. When $`d\sigma _0`$ and $`d\sigma _{\mathrm{HC}}`$ are added together, there is a cancellation between the two singular expressions. Note, however, that this cancellation is not complete since the limits of the $`y`$ integration in the two expressions differ.

After the cancellation, the resulting $`𝒪(\alpha _s)`$ expression for the fragmentation contribution is

$`d\sigma _{frag}^{1+23+h}`$ $`=`$ $`d\sigma _0^{1+23+4^{}}\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]`$
$`\times `$ $`\left\{\stackrel{~}{D}_{h/4^{}}(z,M_f)+\left[{\displaystyle \frac{A_1^{sc}(4^{}4+5)}{ϵ}}+A_0^{sc}(4^{}4+5)\right]D_{h/4}(z,M_f)\right\}dz.`$
The soft collinear factors $`A_i^{sc}`$ result from the mismatch in the $`y`$ integrations in the fragmentation and subtraction pieces, mentioned above. They are given by

$`A_0^{sc}`$ $`=`$ $`A_1^{sc}\mathrm{ln}\left({\displaystyle \frac{s_{12}}{M_f^2}}\right)`$ (61)
$`A_1^{sc}(qqg)`$ $`=`$ $`C_F(2\mathrm{ln}\delta _S+3/2)`$ (62)
$`A_1^{sc}(ggg)`$ $`=`$ $`2N\mathrm{ln}\delta _s+(11N2n_f)/6.`$ (63)
The modified fragmentation function $`\stackrel{~}{D}_{h/c}(z,M_f)`$ is given by

$$\stackrel{~}{D}_{h/c}(z,M_f)=\underset{c^{}}{}_z^{1\delta _s\delta _{c^{}c}}\frac{dy}{y}D_{h/c^{}}(z/y,M_f)\stackrel{~}{P}_{c^{}c}^{\mathrm{frag}},$$
(64)
where

$$\stackrel{~}{P}_{c^{}c}^{\mathrm{frag}}=P_{c^{}c}(y)\mathrm{ln}\left[\frac{y(1y)\delta _cs_{12}}{M_f^2}\right]P_{c^{}c}^{}(y).$$
(65)
$`P(y)`$ and $`P^{}(y)`$ are the $`n=4`$ and $`𝒪(ϵ)`$ pieces, respectively, of the unregulated splitting kernels given in Eqs. (37)–(44). The $`\stackrel{~}{D}`$ functions contain an explicit logarithm of $`\delta _c`$ as well as logarithmic dependences on $`\delta _s`$ which are built up by the integration on $`y`$ when $`c^{}=c`$. In Appendix D it is shown how to make the $`\delta _s`$ dependence explicit, thereby improving the convergence of the Monte Carlo integration.

Comparing with the previous subsection, we see that when going to the fragmentation case, the hard collinear terms, Eqs. (46)–(51), for the fragmenting parton are replaced by a combination of the $`\stackrel{~}{D}`$ function and soft collinear factors $`A_i^{sc}`$. Nevertheless, a careful comparison of the two cases shows that the poles in $`ϵ`$ cancel and the final results for physics observables are independent of the cutoffs. This will be illustrated by several examples to follow.

#### II.2.3 Initial State

The treatment of the initial state collinear singularities is much the same as that for the previous case of final state fragmentation. The collinear singularities are absorbed into the bare parton distribution functions leaving a finite remainder which is written in terms of modified parton distribution functions. In addition, there are accompanying soft collinear factors as in the fragmentation case. However, some of the details are different, so a brief summary of the derivation is given here. In order not to unnecessarily complicate the discussion, only the details for one of the incoming partons will be shown.

Consider a process which involves a parton on leg 2 coming from an incoming hadron $`B`$, so that in lowest order

$$d\sigma _0^{1+B3+4}=G_{2/B}(x)dxd\widehat{\sigma }_0^{1+23+4},$$
(66)
where $`G_{2/B}(x)dx`$ denotes the probability of getting parton 2 from hadron $`B`$ with a momentum fraction between $`x`$ and $`x+dx`$. The hat symbol $`\widehat{}`$ is used here to label a purely partonic subprocess. We are interested in the next-to-leading-order corrections coming from the various possible parton splittings which can occur on leg 2. The hard collinear contribution Eq. (25) is calculated by applying the collinear approximation to the appropriate three-body matrix elements as follows:

$$\overline{}|M_3(1+23+4+5)|^2\overline{}|M_2(1+2^{}3+4)|^2P_{2^{}2}(z,ϵ)g^2\mu _r^{2ϵ}\frac{2}{zt_{25}},$$
(67)
where $`z`$ denotes the fraction of parton 2’s momentum carried by parton $`2^{}`$ with parton 5 taking a fraction $`(1z)`$. Using the approximation $`p_2p_5zp_2`$, the three-body phase space may be written as

$$d\mathrm{\Gamma }_3|_{coll}=\left[\frac{d^{n1}p_3}{2p_3^0(2\pi )^{n1}}\frac{d^{n1}p_4}{2p_4^0(2\pi )^{n1}}(2\pi )^n\delta ^n(p_1+zp_2p_3p_4)\right]\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}.$$
(68)
The square bracketed portion is just two-body phase space evaluated at a squared parton-parton energy of $`zs_{12}`$. The $`p_5`$ dependent part may be rewritten as

$$\frac{d^{n1}p_5}{2p_5^0(2\pi )^{n1}}=\frac{(4\pi )^ϵ}{16\pi ^2\mathrm{\Gamma }(1ϵ)}dzdt_{25}[(1z)t_{25}]^ϵ.$$
(69)
The allowed range for $`t_{25}`$ is given by the collinear condition $`0<t_{25}<\delta _cs_{12}`$. The $`t_{25}`$ integration yields

$$_0^{\delta _cs_{12}}dt_{25}(t_{25})^{1ϵ}=\frac{1}{ϵ}(\delta _cs_{12})^ϵ.$$
(70)
Using these results, the three-body cross section in the hard collinear region may be written as

$`d\sigma _{HC}^{1+B3+4+5}`$ $`=`$ $`G_{2/B}(y)dyd\widehat{\sigma }_0^{1+2^{}3+4}(zs_{12},t_{13},t_{14})\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]`$ (71)
$`\times `$ $`\left({\displaystyle \frac{1}{ϵ}}\right)\delta _c^ϵP_{2^{}2}(z,ϵ)dz(1z)^ϵ\delta (yzx)dx.`$
Note that a factor of $`1/z`$ has been absorbed into the flux factor for the two-body subprocess. The delta function insures that the fraction of hadron $`B^{}s`$ momentum carried by parton $`2^{}`$ into the two-body subprocess is $`x`$ in order to be able to combine this result with the lowest order contribution. The delta function may be used to perform the $`y`$ integration, but one point must first be made. $`s_{12}`$ is related to the square of the overall hadronic squared center-of-mass energy $`S`$ by $`s_{12}=yS`$. On the other hand, in the lowest order subprocess the relation is $`s_{12}=xS`$. It is convenient to rewrite the above expression using this latter definition for $`s_{12}`$. Therefore, after the $`y`$ integration, we obtain

$`d\sigma _{HC}^{1+B3+4+5}`$ $`=`$ $`G_{2/B}(x/z)d\widehat{\sigma }_0^{1+2^{}3+4}(s_{12},t_{13},t_{14})\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]`$ (72)
$`\times `$ $`\left({\displaystyle \frac{1}{ϵ}}\right)\delta _c^ϵP_{2^{}2}(z,ϵ){\displaystyle \frac{dz}{z}}\left[{\displaystyle \frac{(1z)}{z}}\right]^ϵdx.`$
Comparing with the corresponding result for final state fragmentation, we see that a factor of $`[z(1z)]^ϵ`$ has been changed to $`[(1z)/z]^ϵ`$.

In order to factorize the collinear singularity into the parton distribution function, we introduce a scale dependent parton distribution function using the $`\overline{\mathrm{MS}}`$ convention:

$$G_{b/B}(x,\mu _f)=G_{b/B}(x)+\left(\frac{1}{ϵ}\right)\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{\mu _f^2}\right)^ϵ\right]_z^1\frac{dz}{z}P_{bb^{}}(z)G_{b^{}/B}(x/z).$$
(73)
Next, using this definition, we replace $`G_{2/B}(x)`$ in the lowest order expression (66) and combine the result with the hard collinear contribution (72). The resulting $`𝒪(\alpha _s)`$ expression for the initial state collinear contribution is

$`d\sigma _{coll}^{1+B3+4+5}`$ $`=`$ $`d\widehat{\sigma }_0^{1+2^{}3+4}\left[{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\left({\displaystyle \frac{4\pi \mu _r^2}{s_{12}}}\right)^ϵ\right]`$ (74)
$`\times `$ $`\left\{\stackrel{~}{G}_{2^{}/B}(z,\mu _f)+\left[{\displaystyle \frac{A_1^{sc}(22^{}+5)}{ϵ}}+A_0^{sc}(22^{}+5)\right]G_{2/B}(z,\mu _f)\right\}dz.`$
Note that in this expression the soft collinear factors (given in Eqs. (61)–(63)) depend on the initial state factorization scale $`\mu _f`$. The $`\stackrel{~}{G}`$ functions are given by

$$\stackrel{~}{G}_{c/B}(x,\mu _f)=\underset{c^{}}{}_x^{1\delta _s\delta _{cc^{}}}\frac{dy}{y}G_{c^{}/B}(x/y,\mu _f)\stackrel{~}{P}_{cc^{}}(y)$$
(75)
with

$$\stackrel{~}{P}_{ij}(y)=P_{ij}(y)\mathrm{ln}\left(\delta _c\frac{1y}{y}\frac{s_{12}}{\mu _f^2}\right)P_{ij}^{}(y).$$
(76)
The $`n=4`$ and $`𝒪(ϵ)`$ pieces of the unregulated splitting kernels, $`P(y)`$ and $`P^{}(y)`$, are given in Eqs. (37)–(44). An example of a hadron-hadron process will be given in the next section.

As in the final state hadron case, the $`\stackrel{~}{G}`$ functions contain an explicit logarithm of $`\delta _c`$ as well as logarithmic dependences on $`\delta _s`$ which are built up by the integration on $`y`$. In Appendix D it is shown how to make the $`\delta _s`$ dependence explicit, thereby improving the convergence of the Monte Carlo integration.

## III Examples

In this section we provide five illustrative examples applying the method developed in the previous section. The results are shown to be in complete agreement with those available in the literature. We begin by calculating the QCD corrections to electron-positron annihilation into a massive quark pair. The quark mass serves to regulate any would be collinear singularities. There are only final state soft singularities and, hence, only the soft cutoff is required. Next, the QCD corrections to electron-positron annihilation into a massless quark pair are considered. In this case final state soft and collinear singularities are encountered, necessitating the use of both soft and collinear cutoffs. The example of inclusive photon production in hadronic final states of electron-positron annihilation is then presented, illustrating the use of fragmentation functions. Finally, we close the examples section by showing how to calculate the QCD corrections to lepton pair and single particle inclusive production in hadron-hadron collisions. Both examples contain initial state soft and collinear singularities, necessitating the use of scale dependent parton distribution functions. Furthermore, the single particle inclusive cross section calculation also requires the use of scale dependent fragmentation functions.

### III.1 Electron-positron annihilation to massive quark pair

Electron-positron annihilation into a massive quark pair has a particularly simple singularity structure, that of soft singularities in the final state only. It will therefore be used as a first example of the method described in the section above.

Working in the single photon exchange approximation, the leading order Feynman diagram is shown in Fig. 1. Neglecting the electron mass and denoting the quark mass by $`m`$, the leading order cross section

$$d\sigma ^0=\frac{1}{2s_{12}}\overline{}|M_2|^2d\mathrm{\Gamma }_2,$$
(77)
calculated in $`n=42ϵ`$ dimensions is expressed in terms of the (summed and averaged) matrix element squared

$$\overline{}|M_2|^2=2Ne^4Q_q^2\left(\frac{t_{13}^{}{}_{}{}^{2}+t_{23}^{}{}_{}{}^{2}}{s_{12}^2}+\frac{2m^2}{s_{12}}ϵ\right),$$
(78)
and the two-body phase space

$$d\mathrm{\Gamma }_2=\frac{2^{2ϵ}}{16\pi }\left(\frac{4\pi }{s_{12}}\right)^ϵ\beta ^{12ϵ}\frac{1}{\mathrm{\Gamma }(1ϵ)}_0^\pi \mathrm{sin}^{12ϵ}\theta d\theta .$$
(79)
The center-of-mass scattering angle is denoted by $`\theta `$ and $`\beta =\sqrt{14m^2/s_{12}}`$. $`Q_q`$ is the quark charge in units of $`e`$ and $`N=3`$ is the number of colors. When masses are present, it is often convenient to define primed Mandelstam invariants which are the ones defined previously minus some combination of squared masses. In this case, $`s_{ij}^{}s_{ij}m^2`$ and $`t_{ij}^{}t_{ij}m^2`$. Performing the integration over $`\theta `$ we obtain the well known $`ϵ=0`$ result,

$$\sigma ^0=\frac{4\pi \alpha ^2}{3s_{12}}Q_q^2N(1+\frac{2m^2}{s_{12}})\beta ,$$
(80)
with $`\alpha =e^2/4\pi `$.

Because the quark mass regulates any would-be final state collinear singularity, the appropriate decomposition of the two-to-three contribution to the cross section is simply given by Eq. (4). For the $`𝒪(\alpha _s)`$ QCD corrections we therefore need the soft cross section (5), the hard cross section (6), and the virtual corrections.

The real emission diagrams that give $`|M_3|^2`$ are shown in Fig. 2.

If we define the (summed and averaged) squared matrix element as

$$\overline{}|M_3|^2=\frac{1}{4}e^4Q_q^2g^2\frac{1}{s_{12}^2}NC_F\mathrm{\Psi },$$
(81)
then

$`\mathrm{\Psi }`$ $`=`$ $`\mathrm{tr}\left[/p_2\gamma ^\mu /p_1\gamma ^\nu \right]`$ (82)
$`\times `$ $`\{\mathrm{tr}\left[(/p_3+m)\gamma _\mu (/p_4+/p_5m)\gamma ^\sigma (/p_4m)\gamma _\sigma (/p_4/p_5+m)\gamma _\nu \right]/s_{45}^{}{}_{}{}^{2}`$
$`+2\mathrm{tr}\left[(/p_3+m)\gamma ^\sigma (/p_3+/p_5+m)\gamma ^\mu (/p_4m)\gamma _\sigma (/p_4+/p_5m)\gamma _\nu \right]/s_{35}^{}s_{45}^{}`$
$`+\mathrm{tr}\left[(/p_3+m)\gamma ^\sigma (/p_3+/p_5+m)\gamma ^\mu (/p_4+m)\gamma _\nu (/p_3+/p_5+m)\gamma _\sigma \right]/s_{35}^{}{}_{}{}^{2}\},`$
where the strong coupling is denoted by $`g`$ and $`C_F=(N^21)/2N=4/3`$. To obtain the hard contribution, $`\sigma _\mathrm{H}`$, the traces may be evaluated in four space-time dimensions.

There is a soft singularity when the energy of the gluon in Fig. 2 goes to zero. The corresponding soft contribution to the cross section, $`\sigma _\mathrm{S}`$, is given by Eq. (22). In this case, the sum in (22) is taken over the final state quark legs (labeled 3 and 4) and the color connected Born cross sections are related to the leading order cross section by $`d\sigma _{33}^0=d\sigma _{34}^0=d\sigma _{44}^0=C_Fd\sigma ^0`$. We find

$$d\sigma _\mathrm{S}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]C_F\left[\frac{m^2}{(p_3p_5)^2}\frac{m^2}{(p_4p_5)^2}+\frac{s2m^2}{p_3p_5p_4p_5}\right]𝑑S.$$
(83)
The poles need to be integrated over the soft phase space according to Eq. (17) to extract the singularities in dimensional regularization. Define

$`I(s_{35}^{})`$ $`=`$ $`{\displaystyle \frac{1}{s_{35}^{}{}_{}{}^{2}}𝑑S}`$ (84)
$`I(s_{45}^{})`$ $`=`$ $`{\displaystyle \frac{1}{s_{45}^{}{}_{}{}^{2}}𝑑S}`$ (85)
$`I(s_{35}^{}s_{45}^{})`$ $`=`$ $`{\displaystyle \frac{1}{s_{35}^{}s_{45}^{}}𝑑S}.`$ (86)
In terms of the $`p_1p_2`$ center-of-momentum scattering angle $`\theta `$, $`p_3`$ and $`p_4`$ may be written as

$`p_3`$ $`=`$ $`{\displaystyle \frac{\sqrt{s_{12}}}{2}}(1,0,\mathrm{},0,\beta \mathrm{sin}\theta ,\beta \mathrm{cos}\theta )`$
$`p_4`$ $`=`$ $`{\displaystyle \frac{\sqrt{s_{12}}}{2}}(1,0,\mathrm{},0,\beta \mathrm{sin}\theta ,\beta \mathrm{cos}\theta ).`$ (87)
Using these together with Eq. (12) we find

$`s_{35}^{}`$ $`=`$ $`2p_3p_5=\sqrt{s_{12}}E_5(1\beta \mathrm{sin}\theta \mathrm{sin}\theta _1\mathrm{cos}\theta _2\beta \mathrm{cos}\theta \mathrm{cos}\theta _1)`$
$`s_{45}^{}`$ $`=`$ $`2p_4p_5=\sqrt{s_{12}}E_5(1+\beta \mathrm{sin}\theta \mathrm{sin}\theta _1\mathrm{cos}\theta _2+\beta \mathrm{cos}\theta \mathrm{cos}\theta _1).`$ (88)
The gluon energy integrals in Eqs. (84)–(86) may be performed trivially

$$\left(\frac{4}{s_{12}}\right)^ϵ_0^{\delta _s\sqrt{s_{12}}/2}𝑑E_5E_5^{12ϵ}\frac{1}{s_{12}E_5^2}=\frac{1}{s_{12}}\left(\frac{1}{2ϵ}\right)\delta _s^{2ϵ}.$$
(89)
The remaining angular integrals are well know and are tabulated in the Appendix B. The complete results are

$$I(s_{35}^{})=I(s_{45}^{})=\frac{1}{2m^2}\left(\frac{1}{2ϵ}+\mathrm{ln}\delta _s\frac{1}{2\beta }\mathrm{ln}\frac{1+\beta }{1\beta }\right)$$
(90)
$$I(s_{35}^{}s_{45}^{})=\frac{1}{s_{12}\beta }\left(\frac{1}{2ϵ}\mathrm{ln}\frac{1+\beta }{1\beta }\mathrm{Li}_2\frac{2\beta }{1+\beta }\frac{1}{4}\mathrm{ln}^2\frac{1+\beta }{1\beta }+\mathrm{ln}\delta _s\mathrm{ln}\frac{1+\beta }{1\beta }\right).$$
(91)
We may therefore write the final expression for the soft cross section as

$$d\sigma _\mathrm{S}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_1^s}{ϵ}+A_0^s\right),$$
(92)
where

$`A_1^s`$ $`=`$ $`2C_F\left(1{\displaystyle \frac{1+\beta ^2}{2\beta }}\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}\right)`$
$`A_0^s`$ $`=`$ $`4C_F[\mathrm{ln}\delta _s+{\displaystyle \frac{1}{2\beta }}\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}`$ (93)
$``$ $`{\displaystyle \frac{1+\beta ^2}{2\beta }}(\mathrm{Li}_2{\displaystyle \frac{2\beta }{1+\beta }}+{\displaystyle \frac{1}{4}}\mathrm{ln}^2{\displaystyle \frac{1+\beta }{1\beta }}\mathrm{ln}\delta _s\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }})].`$
The virtual contribution is obtained from the one-loop diagrams shown in Fig. 3. In the on-shell renormalization scheme diagrams $`(c)`$ and $`(d)`$ cancel exactly. The vertex correction needed in Fig. 3$`(a)`$ is shown separately in Fig. 4. After performing the loop integrals the result for the vertex valid for $`q^2>4m^2`$ is

$$\mathrm{\Gamma }_\mu =\left(ieQ_q\right)\delta _{ij}\left(g^2C_FC_ϵ\right)\overline{u}(p_2)\left[A\frac{\left(p_1+p_2\right)_\mu }{2m}+B\gamma _\mu \right]u(p_1),$$
(94)
with

$`A`$ $`=`$ $`{\displaystyle \frac{\beta ^21}{\beta }}\mathrm{ln}{\displaystyle \frac{1\beta }{1+\beta }}`$ (95)
$`B`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\left(1+{\displaystyle \frac{1+\beta ^2}{\beta }}\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}\right)+3\beta \mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}`$ (96)
$`+`$ $`{\displaystyle \frac{1+\beta ^2}{\beta }}\left({\displaystyle \frac{1}{2}}\mathrm{ln}^2{\displaystyle \frac{1\beta }{1+\beta }}+2\mathrm{ln}{\displaystyle \frac{1\beta }{1+\beta }}\mathrm{ln}{\displaystyle \frac{2\beta }{1+\beta }}+2\mathrm{Li}_2{\displaystyle \frac{1\beta }{1+\beta }}+{\displaystyle \frac{2\pi ^2}{3}}\right),`$
and

$$C_ϵ=\frac{\pi ^{2ϵ}}{(2\pi )^{42ϵ}}\mathrm{\Gamma }(1+ϵ)\left(\frac{\mu _r^2}{m^2}\right)^ϵ.$$
(97)
The vertex counterterm, $`Z_1`$, implicit in Fig. 3($`b`$) is fixed in the on-shell renormalization scheme by the condition that the renormalized vertex through $`𝒪(g^2)`$ evaluated at zero momentum transfer equals the leading contribution ($`ie\mathrm{\Lambda }^\mu (q=0)=ie\gamma ^\mu `$). This results in

$$Z_1=1g^2C_FC_ϵ\left(\frac{3}{ϵ}+4\right).$$
(98)
The interference of the leading order diagram with the one loop renormalized diagrams yields

$$d\sigma _\mathrm{V}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_1^v}{ϵ}+A_0^v\right)+d\stackrel{~}{\sigma }^v,$$
(99)
with

$`A_1^v`$ $`=`$ $`2C_F\left(1{\displaystyle \frac{1+\beta ^2}{2\beta }}\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}\right)`$
$`A_0^v`$ $`=`$ $`C_F[2(1{\displaystyle \frac{1+\beta ^2}{2\beta }}\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }})\mathrm{ln}\left({\displaystyle \frac{s_{12}}{m^2}}\right)+3\beta \mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}4`$ (100)
$`+`$ $`{\displaystyle \frac{1+\beta ^2}{\beta }}({\displaystyle \frac{1}{2}}\mathrm{ln}^2{\displaystyle \frac{1\beta }{1+\beta }}+2\mathrm{ln}{\displaystyle \frac{1\beta }{1+\beta }}\mathrm{ln}{\displaystyle \frac{2\beta }{1+\beta }}+2\mathrm{Li}_2{\displaystyle \frac{1\beta }{1+\beta }}+{\displaystyle \frac{2}{3}}\pi ^2)],`$
and

$$d\stackrel{~}{\sigma }^v=\frac{8\pi ^2\alpha ^2}{s_{12}}Q_q^2\left(\frac{\alpha _s}{2\pi }\right)\left[4NC_F\frac{\beta ^21}{\beta }\mathrm{ln}\frac{1\beta }{1+\beta }\left(\frac{m^2}{s_{12}}\frac{t_{13}^{}t_{14}^{}}{s_{12}^2}\right)\right]d\mathrm{\Gamma }_2.$$
(101)
Observe that the sum of the soft and ultra-violet renormalized virtual terms is finite, $`A_1^s+A_1^v=0`$, as required kln . We are therefore free to return to 4 dimensions with the finite remainders $`A_0^s+A_0^v`$. The final result for the $`𝒪(\alpha _s)`$ correction consists of two contributions to the cross section: a two-body term $`\sigma ^{(2)}`$ and a three-body term $`\sigma ^{(3)}`$ where

$$\sigma ^{(2)}=\left[d\sigma _0\left(\frac{\alpha _s}{2\pi }\right)\left(A_0^s+A_0^v\right)+d\stackrel{~}{\sigma }^v\right],$$
(102)
and

$$\sigma ^{(3)}=\sigma _H=\frac{1}{2s_{12}}_H\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3.$$
(103)
As a check, these results may be integrated to give a total rate and compared against the known next-to-leading order result zerwas  taken from rn :

$`\sigma _{\mathrm{tot}}^{(1)}`$ $`=`$ $`\sigma _0{\displaystyle \frac{\alpha _s}{4\pi }}C_F[(328\rho ^2)\mathrm{Li}_2(t)+(164\rho ^2)\mathrm{F}_3(t)`$ (104)
$`+`$ $`(2+\rho )\sqrt{1\rho }\mathrm{F}_4(t)+(82\rho ^2)\mathrm{ln}(t)\mathrm{ln}(1+t)`$
$`+`$ $`(12+2\rho +7\rho ^2/4)\mathrm{ln}(t)+(3+9\rho /2)\sqrt{1\rho }],`$
where $`\sigma _0`$ is the leading order cross section for producing a pair of massless quarks given by

$$\sigma _0=\frac{4\pi \alpha ^2}{3s_{12}}NQ_q^2,$$
(105)
and we have also defined

$$t=\frac{1\sqrt{1\rho }}{1+\sqrt{1\rho }},$$
(106)
with $`\rho =4m^2/s_{12}`$ and

$`\mathrm{F}_3(t)`$ $`=`$ $`\mathrm{Li}_2(t)+\mathrm{ln}(t)\mathrm{ln}(1t)`$
$`\mathrm{F}_4(t)`$ $`=`$ $`6\mathrm{ln}(t)8\mathrm{ln}(1t)4\mathrm{ln}(1+t).`$ (107)
The next-to-leading order corrections are shown in Fig. 5 for $`\sqrt{s}=11`$ GeV and $`m=m_b=5`$ GeV. The two- and three-body contributions together with their sum are shown as a function of the soft cutoff $`\delta _s`$. The bottom enlargement shows the sum (open circles) relative to $`\pm 5\%`$ (dotted lines) of the analytical result (solid line) given in Eq. (104). The result quickly converges to the known result.

It is satisfying that the fully inclusive rate from the slicing method agrees with that from Ref. rn . Having made this necessary check, the results may be used to histogram a wide variety observables and to study various physics issues. We refrain from any such studies here, and instead pass to our second example.

### III.2 Electron-positron annihilation to massless quark pair

The process to be studied in this section is similar to that of the last section, but with one key difference: the quarks are considered massless from the beginning. Therefore, in addition to the final state soft singularities there are final state collinear singularities. The leading order cross section

$$d\sigma ^0=\frac{1}{2s_{12}}\overline{}|M_2|^2d\mathrm{\Gamma }_2,$$
(108)
is expressed in terms of the (summed and averaged) matrix element squared

$$\overline{}|M_2|^2=2Ne^4Q_q^2\left(\frac{t_{13}^2+t_{23}^2}{s_{12}^2}ϵ\right),$$
(109)
calculated from Fig. 1, and the two body phase space

$$d\mathrm{\Gamma }_2=\frac{2^{2ϵ}}{16\pi }\left(\frac{4\pi }{s_{12}}\right)^ϵ\frac{1}{\mathrm{\Gamma }(1ϵ)}_0^\pi \mathrm{sin}^{12ϵ}\theta d\theta .$$
(110)
In four space-time dimensions, integration over the phase space produces the result shown previously in Eq. (105).

For the QCD corrections we need the soft cross section (5), the hard-collinear cross section (25), the hard-non-collinear cross section (26), and the virtual contribution. The real emission diagrams are shown in Fig. 2 where the quark lines are to be interpreted as massless. The two-to-three body matrix element squared $`|M_3|^2`$ needed to evaluate the hard-non-collinear cross section Eq. (26) follows directly from Eq. (81) of the previous example by setting $`m=0`$ and evaluating the traces in four space-time dimensions. The soft cross section Eq. (5) may also be obtained from the results of the last example by setting $`m=0`$ in Eq. (83) giving

$$d\sigma _\mathrm{S}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]C_F\left(\frac{s}{p_3p_5p_4p_5}\right)𝑑S.$$
(111)
The pole needs to be integrated over the soft phase space according to Eq. (17). To this end, define

$$I(s_{35}s_{45})=\frac{1}{s_{35}s_{45}}𝑑S.$$
(112)
Using the massless ($`\beta =1`$) form of Eq. (88), the energy integral Eq. (89), and the angular integrals given in Appendix B, the result is

$$I(s_{35}s_{45})=\frac{1}{2s_{12}}\left(\frac{1}{ϵ^2}\frac{2}{ϵ}\mathrm{ln}\delta _s+2\mathrm{ln}^2\delta _s\right).$$
(113)
We may therefore write the final expression for the soft cross section as

$$d\sigma _\mathrm{S}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_2^s}{ϵ^2}+\frac{A_1^s}{ϵ}+A_0^s\right),$$
(114)
with

$`A_2^s`$ $`=`$ $`2C_F`$
$`A_1^s`$ $`=`$ $`4C_F\mathrm{ln}\delta _s`$
$`A_0^s`$ $`=`$ $`4C_F\mathrm{ln}^2\delta _s.`$ (115)
The final state hard collinear cross section was derived in Sec. II.2.1. The relevant splitting is $`qqg`$. From Eq. (45) we have

$$d\sigma _{\mathrm{HC}}^{qqg}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_1^{qqg}}{ϵ}+A_0^{qqg}\right),$$
(116)
with

$`A_1^{qqg}`$ $`=`$ $`C_F\left(3/2+2\mathrm{ln}\delta _s\right)`$
$`A_0^{qqg}`$ $`=`$ $`C_F\left[7/2\pi ^2/3\mathrm{ln}^2\delta _s\mathrm{ln}\delta _c\left(3/2+2\mathrm{ln}\delta _s\right)\right].`$ (117)
The interference of the one-loop diagrams in Fig. 3 with the leading order diagram yields the virtual contribution. In Fig. 3, diagrams $`(b)`$ and $`(d)`$ add to zero via the Ward identity. Diagram $`(c)`$ vanishes for massless quarks. This leaves diagram $`(a)`$, comprised of the vertex shown in Fig. 4 evaluated for massless quarks. The result for the vertex is

$$\mathrm{\Gamma }_\mu =\left(ieQ_q\right)\delta _{ij}\overline{u}(p_2)\gamma _\mu u(p_1)\gamma (q^2),$$
(118)
with

$$\gamma (q^2)=\frac{\alpha _s}{2\pi }C_F\left(\frac{4\pi \mu _r^2}{q^2}\right)^ϵ\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{1}{ϵ^2}+\frac{3}{2ϵ}+4+\frac{\pi ^2}{6}\right).$$
(119)
We may therefore write the final expression for the virtual contribution to the cross section as

$$d\sigma _\mathrm{V}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_2^v}{ϵ^2}+\frac{A_1^v}{ϵ}+A_0^v\right),$$
(120)
with

$`A_2^v`$ $`=`$ $`2C_F`$
$`A_1^v`$ $`=`$ $`3C_F`$
$`A_0^v`$ $`=`$ $`2C_F(4\pi ^2/3).`$ (121)
The full two-body weight is given by the sum $`d\sigma _\mathrm{S}+d\sigma _\mathrm{V}+2d\sigma _{\mathrm{HC}}^{qqg}`$. The factor of two occurs since there are two quark legs, either of which can emit a gluon. At this point we have a finite result since $`A_2^s+A_2^v=0`$ and $`A_1^s+A_1^v+2A_1^{qqg}=0`$ as required kln . The finite two-body weight is given by

$$\sigma ^{(2)}=𝑑\sigma _0\left(\frac{\alpha _s}{2\pi }\right)\left(A_0^s+A_0^v+2A_0^{qqg}\right)$$
(122)
while the three-body contribution is given by

$$\sigma ^{(3)}=\sigma _{H\overline{C}}=\frac{1}{2s_{12}}_{H\overline{C}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3.$$
(123)
A necessary check may be made by integrating these results and comparing with the known analytic answer. The contributions from $`\sigma ^{(2)}`$ (negative) and $`\sigma ^{(3)}`$ (positive) and their sum are shown in Fig. 6 for $`\sqrt{s}=500`$ GeV as a function of the soft cutoff $`\delta _s`$ with the collinear cutoff $`\delta _c=\delta _s/300`$. The known result may be found in sterman , for example, and is given by

$$\sigma _{\mathrm{tot}}^{(1)}=\sigma _0\frac{\alpha _s}{4\pi }3C_F,$$
(124)
where $`\sigma _0`$ is given in Eq. (105). The bottom enlargement shows the sum (open circles) relative to $`\pm 5\%`$ (dotted lines) of the known result (solid line) given in Eq. (124). Very good agreement is found below $`\delta _s2\times 10^3`$.

Before proceeding further, it is instructive to examine some issues related to the cutoff dependence of this technique. As shown in Fig. 6, the answer converges to the known result for $`\delta _s<10^3`$ when $`\delta _c=\delta _s/300`$. We have imposed the requirement $`\delta _c\delta _s`$ which may be understood by examining the nature of the three-body phase space for this case. Neglecting both initial and final state masses, four-momentum conservation yields $`s_{12}=s_{34}+s_{35}+s_{45}`$. The soft region is defined by $`E_5<\delta _s\sqrt{s_{12}}/2`$ which, taken with Eq. (7), can be recast as $`s_{45}<\delta _ss_{12}s_{35}`$. This is shown as the region S in the plot of $`s_{45}`$ versus $`s_{35}`$ in Fig. 7. Two collinear regions defined by the constraints $`s_{35}`$ or $`s_{45}<\delta _cs_{12}`$ are shown as the regions labeled C in Fig. 7. There are two small regions labeled “m” which are properly included in the collinear regions C. However, using a fixed upper limit of $`1\delta _s`$ in calculating the hard collinear contributions (cf. Eq. (35)) these regions are excluded. They are also not included in the hard-non-collinear three-body integrations. With some effort, it is possible to analytically evaluate the required integrals (33) over the m regions. The result (derived in Appendix C) is that occurrences of $`\mathrm{ln}\delta _c\mathrm{ln}\delta _s`$ in Eq. (47) are to be replaced by $`\mathrm{ln}\delta _c\mathrm{ln}\delta _s\mathrm{Li}_2(\delta _c/\delta _s)`$. From the properties of the dilogarithm function we note that the correction term vanishes like $`\delta _c/\delta _s`$ in the limit of small $`\delta _c`$. Accordingly, the contributions from the regions denoted by m in Fig. 7 may be made negligible by requiring $`\delta _c\delta _s`$. Of course, as $`\delta _c`$ and $`\delta _s`$ become smaller the statistical errors on the sum of the two- and three-body weights increase. In practice, one must compromise between the errors induced by larger cutoffs and the statistical errors. For many calculations it has been found that choosing $`\delta _c`$ to be 50 - 100 times smaller than $`\delta _s`$ is sufficient for answers accurate to a few percent. Acceptable ranges for $`\delta _s`$ must be determined on a case by case basis, as illustrated by the examples shown here. Furthermore, the sign of the deviations as $`\delta _s`$ grows differs from process to process.

### III.3 Electron-positron annihilation to photons

In this section we consider an example fragmentation process, inclusive photon production in hadronic final states of electron-positron annihilation, calculated to leading order in the electromagnetic coupling $`\alpha `$. For pedagogical purposes only the radiation of photons from the final state quark or antiquark will be included, i.e., initial state radiation will be neglected. This process is different from the previous two examples in that there are final state collinear singularities only, and they are removed through the factorization procedure.

The diagram for the leading order non-perturbative contribution is shown in Fig. 8. The cross section from Eq. (52) is

$$d\sigma _0^{e^+e^{}\gamma X}=\underset{q}{}dzD_{\gamma /q}(z)d\sigma _0^{e^+e^{}q\overline{q}}.$$
(125)
To simplify notation it is helpful to write the Born-level total cross section for $`e^+e^{}q\overline{q}`$ in terms of that for $`e^+e^{}\mu ^+\mu ^{}`$

$$\sigma _0^{e^+e^{}q\overline{q}}=NQ_q^2\sigma _0^{e^+e^{}\mu ^+\mu ^{}}.$$
(126)
We further denote $`\sigma _0^{e^+e^{}\gamma X}`$ by $`\sigma `$, and $`\sigma _0^{e^+e^{}\mu ^+\mu ^{}}`$ by $`\sigma _{\mu \mu }`$. Taking into account $`D_{\gamma /q}(z)=D_{\gamma /\overline{q}}(z)`$, and letting $`f`$ denote the quark flavor, we arrive at the result for the leading order non-perturbative contribution

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma }{dz}=2N\underset{f=1}{\overset{n_f}{}}Q_f^2D_{\gamma /f}(z).$$
(127)
Additionally, there are two- and three-body pieces that make up the leading order perturbative contribution. The Feynman diagrams are shown in Fig. 9. Because there are only final state collinear singularities present, the relevant decomposition of the two-to-three contribution to the cross section is into collinear and non-collinear terms, $`\sigma =\sigma _\mathrm{C}+\sigma _{\overline{\mathrm{C}}}`$. The collinear term $`\sigma _\mathrm{C}`$ is handled as discussed in Sec. II.2.2. In this case there are no soft singularities so the soft-collinear terms $`A^{sc}`$ are not present. The two-body piece follows from Eq. (II.2.2)

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma ^{(2)}}{dz}=2N\frac{\alpha }{2\pi }\underset{f=1}{\overset{n_f}{}}Q_q^4\stackrel{~}{D}_{\gamma /f}(z,M_f),$$
(128)
once the replacement $`\alpha _s\alpha `$ is made. $`\stackrel{~}{D}_{\gamma /f}(z,M_f)`$ is as given in Eq. (64) and may be expanded using $`D_{\gamma /\gamma }(x)=\delta (1x)+𝒪(\alpha ^2)`$ and $`D_{\gamma /i}(x)=𝒪(\alpha )`$ for $`i=g,q`$. The leading term in $`\alpha `$ is therefore

$$\stackrel{~}{D}_{\gamma /q}(z,M_f)=_z^1\frac{dy}{y}\delta (1z/y)\stackrel{~}{P}_{\gamma q}^{\mathrm{frag}}(y)=\stackrel{~}{P}_{\gamma q}^{\mathrm{frag}}(z),$$
(129)
where $`\stackrel{~}{P}^{\mathrm{frag}}`$ is as given in Eq. (65) with

$$P_{\gamma q}(z,ϵ)=\frac{1+(1z)^2ϵz^2}{z}.$$
(130)
The final result for the two-body piece of the leading order perturbative contribution is

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma ^{(2)}}{dz}=2N\frac{\alpha }{2\pi }\underset{f=1}{\overset{n_f}{}}Q_f^4\left\{\frac{1+(1z)^2}{z}\mathrm{ln}\left[\frac{z(1z)\delta _cs}{M_f^2}\right]+z\right\}.$$
(131)
The complementary non-collinear piece $`\sigma _{\overline{\mathrm{C}}}`$ follows from the matrix element represented in Fig. 9

$$\sigma ^{(3)}=\frac{1}{2s}_{\overline{\mathrm{C}}}\overline{}|M_3|^2𝑑\mathrm{\Gamma }_3=\frac{4e^6}{3s^2}N\underset{f=1}{\overset{n_f}{}}Q_f^4_{\overline{\mathrm{C}}}𝑑\mathrm{\Gamma }_3\left(\frac{2s_{34}s}{s_{35}s_{45}}+\frac{s_{35}}{s_{45}}+\frac{s_{45}}{s_{35}}\right).$$
(132)
We could now study the cutoff dependence numerically as in the previous examples. However, the integration over phase space may be performed analytically and rather straightforwardly so we take this route to demonstrating the $`\delta _c`$ independence of the full result. To this end, consider, in the virtual photon rest frame, the three-body phase space

$$d\mathrm{\Gamma }_3=\frac{1}{8}\frac{1}{(2\pi )^5}\frac{1}{E_3}\frac{d^3p_4}{E_4}\frac{d^3p_5}{E_5}\delta (\sqrt{s}E_3E_4E_5).$$
(133)
Let $`q^\mu `$ denote the virtual photon four-momentum. Taking $`\stackrel{}{p}_5`$ along the $`z`$ axis and defining $`x_4=2E_4/\sqrt{s}`$, $`z=2E_5/\sqrt{s}`$ we may write the four-momenta as

$`q`$ $`=`$ $`\sqrt{s}(1,0,0,0)`$ (134)
$`p_5`$ $`=`$ $`z{\displaystyle \frac{\sqrt{s}}{2}}(1,0,0,1)`$ (135)
$`p_4`$ $`=`$ $`x_4{\displaystyle \frac{\sqrt{s}}{2}}(1,\mathrm{sin}\theta ,0,\mathrm{cos}\theta ).`$ (136)
Momentum conservation $`\stackrel{}{q}=\stackrel{}{p}_3+\stackrel{}{p}_4+\stackrel{}{p}_5`$ gives

$$p_3=(E_3,x_4\frac{\sqrt{s}}{2}\mathrm{sin}\theta ,0,x_4\frac{\sqrt{s}}{2}\mathrm{cos}\theta z\frac{\sqrt{s}}{2}).$$
(137)
The mass-shell condition $`p_3^2=0`$ can be used to fix $`E_3`$ as

$$E_3=\frac{\sqrt{s}}{2}\sqrt{x_4^2+z^2+2x_4z\mathrm{cos}\theta }.$$
(138)
Writing $`d^3p_4=2\pi d\mathrm{cos}\theta E_4^2dE_4`$ and $`d^3p_5=4\pi E_5^2dE_5`$ the phase space delta function may be used to perform the $`\mathrm{cos}\theta `$ integral. The invariants $`s_{34}`$, $`s_{35}`$, and $`s_{45}`$ may then be written in terms of $`x_4`$ and $`z`$

$`s_{35}`$ $`=`$ $`s(1x_4)`$
$`s_{45}`$ $`=`$ $`s(z+x_41)`$
$`s_{34}`$ $`=`$ $`s(1z).`$ (139)
The three-body piece is now

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma ^{(3)}}{dz}=N\frac{\alpha }{2\pi }\underset{f=1}{\overset{n_f}{}}Q_f^4_{\overline{\mathrm{C}}}𝑑x_4\left[\frac{2(1z)}{(1x_4)(z+x_41)}+\frac{1x_4}{z+x_41}+\frac{z+x_41}{1x_4}\right].$$
(140)
This is to be integrated over the non-collinear $`\overline{\mathrm{C}}`$ region defined by $`s_{45}>\delta _cs`$ and $`s_{35}>\delta _cs`$ which is equivalent to $`1z+\delta _cx_41\delta _c`$. The integral may easily be performed. Dropping terms of $`𝒪(\delta _c)`$, the final result for the three-body piece of the leading order perturbative contribution is

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma ^{(3)}}{dz}=2N\frac{\alpha }{2\pi }\underset{f=1}{\overset{n_f}{}}Q_f^4\left[\frac{1+(1z)^2}{z}\mathrm{ln}\left(\frac{z}{\delta _c}\right)z\right].$$
(141)
Adding the non-perturbative (127), perturbative two-body (131), and perturbative three-body (141) contributions we obtain the well known result photon1 ; photon2

$$\frac{1}{\sigma _{\mu \mu }}\frac{d\sigma }{dz}=2N\underset{f=1}{\overset{n_f}{}}Q_f^2\left\{D_{\gamma /f}(z)+\frac{\alpha }{2\pi }Q_f^2\frac{1+(1z)^2}{z}\mathrm{ln}\left[\frac{z^2(1z)s}{M_f^2}\right]\right\},$$
(142)
which is independent of $`\delta _c`$.

### III.4 Drell-Yan

Our next example is that of the QCD corrections to lepton pair production in hadron-hadron collisions which illustrates the method for handling initial state collinear singularities developed in Sec. II.2.3.

The leading order contribution mediated by a virtual photon is shown in Fig. 10. The leading order partonic cross section

$$d\widehat{\sigma }^0=\frac{1}{2s_{12}}\overline{}|M_2|^2d\mathrm{\Gamma }_2,$$
(143)
is expressed in terms of the (summed and averaged) matrix element squared calculated in $`n=42ϵ`$ dimensions

$$\overline{}|M_2|^2=e^4Q_f^2\frac{2}{N}\left(\frac{t_{13}^2+t_{23}^2}{s^2}ϵ\right),$$
(144)
and the two-body phase space

$$d\mathrm{\Gamma }_2=\frac{2^{2ϵ}}{16\pi }\left(\frac{4\pi }{s_{12}}\right)^ϵ\frac{1}{\mathrm{\Gamma }(1ϵ)}_0^\pi \mathrm{sin}^{12ϵ}\theta d\theta .$$
(145)
For the QCD corrections there are five pieces to consider: the finite hard-non-collinear partonic cross section Eq. (26), the mass factorized hard-collinear cross section Eq. (74), which consists of two pieces, the soft part of the initial state factorization counterterms and the mass factorization residuals (the $`\stackrel{~}{G}`$ functions), the soft cross section Eq. (22), and the virtual corrections.

Shown in Fig. 11 are the real emission diagrams that give $`|M_3|^2`$. Defining the (summed and averaged) matrix element squared as

$$\overline{}|M_3|^2=\frac{1}{4}e^4Q_f^2g^2\frac{1}{s_{34}^2}\frac{C_F}{N}\mathrm{\Psi }.$$
(146)
we find

$`\mathrm{\Psi }`$ $`=`$ $`\mathrm{tr}\left[/p_3\gamma ^\mu /p_4\gamma ^\nu \right]`$ (147)
$`\times `$ $`\{\mathrm{tr}\left[/p_2\gamma _\mu (/p_1/p_5)\gamma ^\sigma /p_1\gamma _\sigma (/p_1/p_5)\gamma _\nu \right]/t_{15}^2`$
$`+2\mathrm{tr}\left[/p_2\gamma ^\sigma (/p_2/p_5)\gamma ^\mu /p_1\gamma _\sigma (/p_1/p_5)\gamma _\nu \right]/t_{15}t_{25}`$
$`\mathrm{tr}\left[/p_2\gamma ^\sigma (/p_2/p_5)\gamma ^\mu /p_1\gamma _\nu (/p_2/p_5)\gamma _\sigma \right]/t_{25}^2\}.`$
The hard-non-collinear partonic cross section is obtained by evaluating the traces in four space-time dimensions.

There is a soft singularity when the gluon’s energy goes to zero in the real emission diagrams. This contributes to the soft cross section presented in Eq. (22). The sum runs over the initial state quark lines (labeled by 1 and 2). The color connected Born cross section $`d\sigma _{12}^0=C_Fd\sigma ^0`$. We find

$$d\sigma _\mathrm{S}=d\widehat{\sigma }^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]C_F\left(\frac{s}{p_1p_5p_2p_5}\right)𝑑S.$$
(148)
The integration of the pole over the soft phase space measure (17) is written in terms of

$$I(t_{15}t_{25})=\frac{1}{t_{15}t_{25}}𝑑S.$$
(149)
In the $`p_1p_2`$ center-of-momentum system we take

$`p_1`$ $`=`$ $`{\displaystyle \frac{\sqrt{s_{12}}}{2}}(1,0,\mathrm{},0,1)`$
$`p_2`$ $`=`$ $`{\displaystyle \frac{\sqrt{s_{12}}}{2}}(1,0,\mathrm{},0,1).`$ (150)
Using these together with Eq. (12) we find

$`t_{15}`$ $`=`$ $`2p_1p_5=\sqrt{s_{12}}E_5(1\mathrm{cos}\theta _1)`$
$`t_{25}`$ $`=`$ $`2p_2p_5=\sqrt{s_{12}}E_5(1+\mathrm{cos}\theta _1).`$ (151)
Using the energy integral Eq. (89) and the angular integrals given in Appendix B, we find

$$I(t_{15}t_{25})=\frac{1}{2s_{12}}\left(\frac{1}{ϵ^2}\frac{2}{ϵ}\mathrm{ln}\delta _s+2\mathrm{ln}^2\delta _s\right).$$
(152)
We may therefore write the final expression for the soft cross section as

$$d\sigma _\mathrm{S}=d\widehat{\sigma }^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_2^s}{ϵ^2}+\frac{A_1^s}{ϵ}+A_0^s\right),$$
(153)
with

$`A_2^s`$ $`=`$ $`2C_F`$
$`A_1^s`$ $`=`$ $`4C_F\mathrm{ln}\delta _s`$
$`A_0^s`$ $`=`$ $`4C_F\mathrm{ln}^2\delta _s.`$ (154)
Because the quarks are massless, there is a collinear singularity when the gluon becomes collinear to either of the initial state quark lines. This singularity is removed through the factorization program described in Sec. II.2.3. The soft-collinear pieces of the initial state factorization counterterms are given in Eq. (74) as

$$d\sigma ^{\mathrm{coll}}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_1^{\mathrm{sc}}}{ϵ}+A_0^{\mathrm{sc}}\right),$$
(155)
with

$`A_1^{\mathrm{sc}}`$ $`=`$ $`C_F\left(2\mathrm{ln}\delta _s+3/2\right)`$
$`A_0^{\mathrm{sc}}`$ $`=`$ $`C_F\left(2\mathrm{ln}\delta _s+3/2\right)\mathrm{ln}{\displaystyle \frac{s_{12}}{\mu _f^2}}.`$ (156)
The one-loop virtual diagrams are shown in Fig. 12. As in the case for electron-positron annihilation to a massless quark pair, diagrams $`(b)`$ and $`(d)`$ add to zero via the Ward identity, and diagram $`(c)`$ vanishes for massless quarks. This leaves diagram $`(a)`$ for which the vertex shown in Fig. 4 is needed, for massless quarks. The result is given in Eq. (118). The final expression for the virtual cross section is

$$d\sigma _\mathrm{V}=d\sigma ^0\left[\frac{\alpha _s}{2\pi }\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}\left(\frac{4\pi \mu _r^2}{s_{12}}\right)^ϵ\right]\left(\frac{A_2^v}{ϵ^2}+\frac{A_1^v}{ϵ}+A_0^v\right),$$
(157)
with

$`A_2^v`$ $`=`$ $`2C_F`$
$`A_1^v`$ $`=`$ $`3C_F`$
$`A_0^v`$ $`=`$ $`2C_F(4\pi ^2/3).`$ (158)
At this point we pause to note that the two-body weight is finite: $`A_2^s+A_2^v=0`$ and $`A_1^s+A_1^v+2A_1^{\mathrm{sc}}=0`$. The factor of two occurs since there are two quark legs, either of which can emit a gluon.

In addition to the $`q\overline{q}`$ initiated processes, there are also $`qg`$ initiated processes at this order of perturbation theory, as shown in Fig. 13. The singularities are initial state collinear only in origin and arise from the $`P_{qg}`$ splitting in diagram $`(b)`$. They are removed by factorization. As the $`P_{qg}(z)`$ kernel is finite for $`z=1`$ there are no soft singularities. This implies the $`A^{\mathrm{sc}}`$ terms in Eq. (74) are not present; only the finite $`\stackrel{~}{G}`$ terms remain.

The final finite two-body cross section is given by the sum of the residual $`\stackrel{~}{G}`$ terms from both the $`q\overline{q}`$ and $`qg`$ initiated processes and the finite two body weights from the $`q\overline{q}`$ process, $`A_0^s+A_0^v+2A_0^{\mathrm{sc}}`$. The result, summed over all parton flavors, is

$`\sigma ^{(2)}`$ $`=`$ $`\left({\displaystyle \frac{\alpha _s}{2\pi }}\right){\displaystyle \underset{f}{}}{\displaystyle }dx_Adx_Bd\widehat{\sigma }_0[G_{f/A}(x_A,\mu _f)G_{\overline{f}/B}(x_B,\mu _f)(A_0^s+A_0^v+2A_0^{\mathrm{sc}})`$ (159)
$`+`$ $`G_{f/A}(x_A,\mu _f)\stackrel{~}{G}_{\overline{f}/B}(x_B,\mu _f)+\stackrel{~}{G}_{f/A}(x_A,\mu _f)G_{\overline{f}/B}(x_B,\mu _f)+(x_Ax_B)].`$
The $`\stackrel{~}{G}`$ functions are given in Eq. (75). The three-body contribution is given by

$$\sigma ^{(3)}=\underset{i,j=q,g}{}𝑑x_A𝑑x_BG_{i/A}(x_A,\mu _f)G_{j/B}(x_B,\mu _f)𝑑\widehat{\sigma }_{ij},$$
(160)
with the hard-non-collinear partonic cross section given by

$$d\widehat{\sigma }_{ij}=\frac{1}{2s_{12}}_{H\overline{C}}\overline{}|M_3^{(ij)}|^2𝑑\mathrm{\Gamma }_3.$$
(161)
Physical predictions follow from the sum $`\sigma ^{(2)}+\sigma ^{(3)}`$ which is cutoff independent for sufficiently small cutoffs. The results may be integrated to obtain the total rate for $`Q^2>Q_{\mathrm{min}}^2`$ and checked against the known $`𝒪(\alpha _s)`$ corrections

$`\sigma ={\displaystyle _{Q_{\mathrm{min}}^2}^S}dQ^2{\displaystyle _{Q^2/S}^1}dx_A{\displaystyle _{Q^2/Sx_A}^1}dx_B\{{\displaystyle \underset{ij=q,\overline{q}}{}}G_{i/A}(x_A,\mu _f)G_{j/B}(x_B,\mu _f){\displaystyle \frac{d\widehat{\sigma }_{qq}}{dQ^2}}`$ (162)
$`+{\displaystyle \underset{i=q,\overline{q}}{}}[G_{i/A}(x_A,\mu _f)G_{g/B}(x_B,\mu _f)+G_{g/A}(x_A,\mu _f)G_{i/B}(x_B,\mu _f)]{\displaystyle \frac{d\widehat{\sigma }_{qg}}{dQ^2}}\},`$
where $`Q^2`$ is the square of the lepton pair invariant mass and $`S`$ is the hadron-hadron center of mass energy squared which is related to $`s_{12}`$, the parton-parton center of mass energy squared, via $`s_{12}=x_Ax_BS`$. Defining $`z=Q^2/s_{12}`$, the $`𝒪(\alpha _s)`$ hard scattering partonic subprocess cross sections are given by aem

$`{\displaystyle \frac{d\widehat{\sigma }_{qq}}{dQ^2}}`$ $`=`$ $`\sigma _0{\displaystyle \frac{\alpha _s}{2\pi }}C_F\{4(1+z^2)\left[{\displaystyle \frac{\mathrm{ln}(1z)}{1z}}\right]_+2{\displaystyle \frac{(1+z^2)}{1z}}\mathrm{ln}z`$ (163)
$`+({\displaystyle \frac{2\pi ^2}{3}}8)\delta (1z)+{\displaystyle \frac{3}{2}}P_{qq}^+(z)\mathrm{ln}{\displaystyle \frac{Q^2}{\mu _f^2}}\},`$
and

$`{\displaystyle \frac{d\widehat{\sigma }_{qg}}{dQ^2}}`$ $`=`$ $`\sigma _0{\displaystyle \frac{\alpha _s}{2\pi }}{\displaystyle \frac{1}{2}}\left\{{\displaystyle \frac{3}{2}}+z{\displaystyle \frac{3}{2}}z^2+2P_{qg}^+(z)\left[\mathrm{ln}{\displaystyle \frac{(1z)^2}{z}}1+\mathrm{ln}{\displaystyle \frac{Q^2}{\mu _f^2}}\right]\right\},`$ (164)
where $`\sigma _0`$

$$\sigma _0=\frac{4\pi \alpha ^2Q_q^2}{3Ns_{12}Q^2}.$$
(165)
$`P_{qq}^+(z)`$ and $`P_{qg}^+(z)`$ are the splitting kernels given in Eqs. (55) and (58).

We show numerical results on the cutoff (in)dependence for proton-proton collisions at $`\sqrt{S}=28.28`$ GeV with $`Q_{min}=10`$ GeV. Hard scales are set to the lepton pair mass $`\mu _f=\mu _r=Q`$ and the number of flavors taken to be $`n_f=3`$.

Shown in Fig. 14 is the next-to-leading order quark-quark contribution to the Drell-Yan cross section. The two- and three-body contributions to the cross section (negative and positive, respectively), and their sum are shown as a function of the soft cutoff $`\delta _s`$. The collinear cutoff $`\delta _c=\delta _s/50`$. The bottom portion of the figure shows the sum (open circles) relative to $`\pm 5\%`$ (dotted lines) of the analytic result (solid line) given in Eq. (163).

Finally, we show the next-to-leading order quark-gluon contribution to the Drell-Yan cross section in Fig. 15. The two- and three-body contributions and their sum are shown as a function of the collinear cutoff $`\delta _c`$. The bottom enlargement shows the sum (open circles) relative to $`\pm 5\%`$ (dotted lines) of the analytic result (solid line) given in Eq. (164).

In both cases, nice agreement is seen with the known analytic result, providing a cross check on the use of the two cutoff phase space slicing method.

### III.5 Single Particle Inclusive Cross Section

Our final example is that of the single particle inclusive cross section in hadron-hadron collisions. The input needed for this calculation includes the squared matrix elements for the $`23`$ subprocesses es  and the results for the $`𝒪(\alpha _s^3)`$ one-loop contributions to the $`22`$ subprocesses ks ; es . For the purpose of this example, the notation of ks  will be used, since much of the input needed can be found in the appendices of that paper. The partons are labelled as $`A+B1+2`$ and $`A+B1+2+3`$ for the $`22`$ and $`23`$ subprocesses, respectively. A flavor label $`a_A`$ is used to denote the flavor of parton $`A`$, and similarly for the other partons.

The lowest-order contribution to the inclusive cross section for producing a hadron $`h`$ in a collision of hadrons of types $`A`$ and $`B`$ can be written as

$`d\sigma ^B`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}G_{a_A/A}(x_A)G_{a_B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (166)
$`\times `$ $`{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})d\mathrm{\Gamma }_2`$
where $`\stackrel{}{a}=\{a_A,a_B,a_1,a_2\}`$ and $`\stackrel{}{p}=\{p_A^\mu ,p_B^\mu ,p_1^\mu ,p_2^\mu \}`$ denote the sets of flavor indices and parton four-vectors, respectively. The factors appearing in the spin/color averaging are given by

$$w(a)=\{\begin{array}{cc}2(1ϵ)V\hfill & \text{a=gluon}\hfill \\ 2N\hfill & \text{a=quark or antiquark}\hfill \end{array}$$

with $`N=3`$ and $`V=N^21.`$ The factor $`d\mathrm{\Gamma }_2`$ is the differential two-body phase space element from Eq. (11). Eq. (166) gives the contribution where parton 1 fragments into the hadron $`h`$. Care must be taken to explicitly include in the sum over $`\stackrel{}{a}`$ those terms corresponding to the case where parton 2 fragments into $`h`$. For compactness, these terms will not be explicitly written. The squared matrix elements for the various subprocesses, denoted by $`\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})`$, may be found in Ref. ks .

Next, consider the one-loop virtual corrections to the $`22`$ subprocesses. These take the form

$`d\sigma ^v`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}G_{a_A/A}(x_A)G_{a_B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (167)
$`\times `$ $`{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}\left[{\displaystyle \frac{\alpha _s}{2\pi }}\left({\displaystyle \frac{4\pi \mu _R^2}{2p_Ap_B}}\right)^ϵ{\displaystyle \frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}}\right]\psi ^{(6)}(\stackrel{}{a},\stackrel{}{p})d\mathrm{\Gamma }_2`$
where

$`\psi ^{(6)}(\stackrel{}{a},\stackrel{}{p})`$ $`=`$ $`\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})\left[{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \underset{n}{}}C(a_n){\displaystyle \frac{1}{ϵ}}{\displaystyle \underset{n}{}}\gamma (a_n)\right]`$ (168)
$`+`$ $`{\displaystyle \frac{1}{2ϵ}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{m,n}{mn}}{}}\mathrm{ln}\left({\displaystyle \frac{p_mp_n}{p_Ap_B}}\right)\psi _{m,n}^{(4,c)}(\stackrel{}{a},\stackrel{}{p})`$
$``$ $`{\displaystyle \frac{\pi ^2}{6}}{\displaystyle \underset{n}{}}\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})+\psi _{NS}^{(6)}(\stackrel{}{a},\stackrel{}{p})+𝒪(ϵ).`$
This expression for $`\psi ^{(6)}`$ differs slightly from Eq. (35) in Ref. ks  because we have chosen to extract a different $`ϵ`$ dependent overall factor: a factor of $`\mathrm{\Gamma }(1+ϵ)\mathrm{\Gamma }(1ϵ)1+ϵ^2\frac{\pi ^2}{6}`$ has been absorbed into the above expression for $`\psi ^{(6)}`$. Furthermore, the arbitrary scale $`Q_{ES}^2`$ used in Ref. ks  has, in order to match the conventions used elsewhere in the present work, been chosen to be $`2p_Ap_B`$. The expressions for the functions $`\psi _{m,n}^{(4,c)}`$ and $`\psi _{NS}^{(6)}`$ may be found in Appendix B of Ref. ks . The quantities $`C(a_n)`$ and $`\gamma (a_n)`$ are given by

$$C(a)=\{\begin{array}{cc}N=3\hfill & \text{a=gluon}\hfill \\ C_F=\frac{4}{3}\hfill & \text{a=quark or antiquark}\hfill \end{array}$$

and

$$\gamma (a)=\{\begin{array}{cc}(11N2n_f)/6\hfill & \text{a=gluon}\hfill \\ 3C_F/2\hfill & \text{a=quark or antiquark}\hfill \end{array}$$

It will be convenient for subsequent expressions to adopt the following notation:

$$=\left(\frac{4\pi \mu _R^2}{2p_Ap_B}\right)^ϵ\frac{\mathrm{\Gamma }(1ϵ)}{\mathrm{\Gamma }(12ϵ)}.$$
(169)
The one loop virtual contributions can now be written as

$`d\sigma ^v`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}G_{a_A/A}(x_A)G_{a_B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (170)
$`\times `$ $`{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}{\displaystyle \frac{\alpha _s}{2\pi }}\left({\displaystyle \frac{A_2^v}{ϵ^2}}+{\displaystyle \frac{A_1^v}{ϵ}}+A_0^v\right)d\mathrm{\Gamma }_2`$
where

$`A_2^v`$ $`=`$ $`{\displaystyle \underset{n}{}}C(a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})`$ (171)
$`A_1^v`$ $`=`$ $`{\displaystyle \underset{n}{}}\gamma (a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{m,n}{mn}}{}}\mathrm{ln}\left({\displaystyle \frac{p_mp_n}{p_Ap_B}}\right)\psi _{m,n}^{(4,c)}(\stackrel{}{a},\stackrel{}{p})`$ (172)
$`A_0^v`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{6}}{\displaystyle \underset{n}{}}C(a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})+\psi _{NS}^{(6)}(\stackrel{}{a},\stackrel{}{p}).`$ (173)
Next, the contributions from the $`23`$ subprocesses in the limit where one of the final state gluons becomes soft are needed. The contributions of the $`23`$ subprocesses may be written as

$`d\sigma ^{23}`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \frac{(4\pi \alpha _s)^3}{w(a_A)w(a_B)}}{\displaystyle \underset{a_A,a_B,a_1,a_2,a_3}{}}G_{a_A/A}(x_A)G_{a/B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (174)
$`\times `$ $`\mathrm{\Psi }(a_A,a_B,a_1,a_2,a_3,p_A^\mu ,p_B^\mu ,p_1^\mu ,p_2^\mu ,p_3^\mu )d\mathrm{\Gamma }_3.`$
The expressions for the $`23`$ squared matrix elements appearing in Eq. (174) may be found in Ref. es . As noted earlier for the two-body contributions, one must include in the sum all possible parton to hadron fragmentations.

Consider the case where the soft gluon is parton 3. In this limit, the function $`\mathrm{\Psi }`$ may be expanded as:

$$\mathrm{\Psi }(a_A,a_B,a_1,a_2,a_3,p_A^\mu ,p_B^\mu ,p_1^\mu ,p_2^\mu ,p_3^\mu )\underset{\genfrac{}{}{0pt}{}{m,n}{m<n}}{}\delta _{a_3,g}\frac{p_mp_n}{p_mp_3p_np_3}\psi _{m,n}^{(4,c)}(a_A,a_B,a_1,a_2,p_A^\mu ,p_B^\mu ,p_1^\mu ,p_2^\mu ).$$
(175)
Next, one must integrate over the soft region of phase space defined by $`E_3<\delta _s\sqrt{2p_Ap_B}/2`$. This is easily done using the integrals given in the appendix. The resulting soft contribution may be written as

$`d\sigma ^s`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}G_{a_A/A}(x_A)G_{a_B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (176)
$`\times `$ $`{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}{\displaystyle \frac{\alpha _s}{2\pi }}\left({\displaystyle \frac{A_2^s}{ϵ^2}}+{\displaystyle \frac{A_1^s}{ϵ}}+A_0^s\right)d\mathrm{\Gamma }_2`$
where

$`A_2^s`$ $`=`$ $`{\displaystyle \underset{n}{}}C(a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})`$ (177)
$`A_1^s`$ $`=`$ $`2\mathrm{ln}\delta _s{\displaystyle \underset{n}{}}C(a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p}){\displaystyle \frac{1}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{m,n}{mn}}{}}\mathrm{ln}\left({\displaystyle \frac{p_mp_n}{p_Ap_B}}\right)\psi _{m,n}^{(4,c)}(\stackrel{}{a},\stackrel{}{p})`$ (178)
$`A_0^s`$ $`=`$ $`2\mathrm{ln}^2\delta _s{\displaystyle \underset{n}{}}C(a_n)\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})`$ (179)
$`+`$ $`\left(\psi _{A,1}^{(4,c)}+\psi _{B,2}^{(4,c)}\right)\left[{\displaystyle \frac{1}{2}}\mathrm{ln}^2\left({\displaystyle \frac{p_1p_3}{p_Ap_B}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{p_2p_3}{p_Ap_B}}\right)+2\mathrm{ln}\delta _s\mathrm{ln}\left({\displaystyle \frac{p_1p_3}{p_Ap_B}}\right)\right]`$
$`+`$ $`\left(\psi _{A,2}^{(4,c)}+\psi _{B,1}^{(4,c)}\right)\left[{\displaystyle \frac{1}{2}}\mathrm{ln}^2\left({\displaystyle \frac{p_2p_3}{p_Ap_B}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{p_1p_3}{p_Ap_B}}\right)+2\mathrm{ln}\delta _s\mathrm{ln}\left({\displaystyle \frac{p_2p_3}{p_Ap_B}}\right)\right].`$
After the collinear singularities associated with the two parton distribution functions and the fragmentation function have been factorized and absorbed into the corresponding bare functions, there will be soft-collinear terms left over due to the mismatch between the integration limits of the collinear singularity terms and the factorization counterterms. In addition, there can be collinear singularities associated with the non-fragmenting parton in the final state, corresponding to gluon emission or $`q\overline{q}`$ production. Collecting together both types of collinear terms, the result can be written as follows:

$`d\sigma ^{coll}`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}G_{a_A/A}(x_A)G_{a_B/B}(x_B)D_{h/a_1}(z_1)dx_Adx_Bdz_1`$ (180)
$`\times `$ $`{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}{\displaystyle \frac{\alpha _s}{2\pi }}\left({\displaystyle \frac{A_1^{coll}}{ϵ}}+A_0^{coll}\right)d\mathrm{\Gamma }_2`$
where

$`A_1^{coll}`$ $`=`$ $`{\displaystyle \underset{n}{}}[2\mathrm{ln}\delta _sC(a_n)+\gamma (a_n)]`$ (181)
$`A_0^{coll}`$ $`=`$ $`{\displaystyle \underset{A,B}{}}[2\mathrm{ln}\delta _sC(a_n)+\gamma (a_n)]\mathrm{ln}\left({\displaystyle \frac{2p_Ap_B}{\mu _f^2}}\right)+[2\mathrm{ln}\delta _sC(a_1)+\gamma (a_1)]\mathrm{ln}\left({\displaystyle \frac{2p_Ap_B}{M_f^2}}\right)`$ (182)
$`+`$ $`\gamma ^{}(a_2).`$
Here $`\mu _f\text{ and }M_f`$ are the initial and final state factorization scales. The function $`\gamma ^{}(a)`$ is given in terms of the hard collinear factors of Eqs. (46) – (51) as

$$\gamma ^{}(a)=\{\begin{array}{cc}A_0^{qqg}\hfill & \text{a= quark or antiquark}\hfill \\ A_0^{ggg}+A_0^{gq\overline{q}}\hfill & \text{a = gluon}\hfill \end{array}$$
(183)
After the mass factorization has been performed, the bare parton distribution functions and fragmentation functions have been replaced by scale dependent $`\overline{\text{MS}}`$ functions. In addition, there are finite remainders involving the $`\stackrel{~}{G}\text{ and }\stackrel{~}{D}`$ functions:

$`d\stackrel{~}{\sigma }`$ $`=`$ $`{\displaystyle \frac{1}{2x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}{\displaystyle \frac{\alpha _s}{2\pi }}dx_Adx_Bdz_1\psi ^{(4)}(\stackrel{}{a},\stackrel{}{p})d\mathrm{\Gamma }_2`$ (184)
$`\times `$ $`[\stackrel{~}{G}_{a_A/A}(x_A,\mu _f^2)G_{a_B/B}(x_B,\mu _f^2)D_{h/a_1}(z_1,M_f^2)`$
$`+`$ $`G_{a_A/A}(x_A,\mu _f^2)\stackrel{~}{G}_{a_B/B}(x_B,\mu _f^2)D_{h/a_1}(z_1,M_f^2)`$
$`+`$ $`G_{a_A/A}(x_A,\mu _f^2)G_{a_B/B}(x_B,\mu _f^2)\stackrel{~}{D}_{h/a_1}(z_1,M_f^2)].`$
At this point, all of the singular terms have been isolated as poles in $`ϵ`$ or have been factorized and absorbed into the bare parton distribution and fragmentation functions. The $`ϵ`$ dependent pole terms all cancel amongst each other:

$`A_2^v+A_2^s=0`$ (185)
$`A_1^v+A_1^s+A_1^{coll}=0.`$ (186)
The finite two-body contribution is given by

$`d\sigma ^{22}`$ $`=`$ $`d\sigma ^B+d\stackrel{~}{\sigma }`$ (187)
$`+`$ $`{\displaystyle \frac{1}{x_Ax_Bs}}{\displaystyle \underset{a_A,a_B,a_1,a_2}{}}{\displaystyle \frac{(4\pi \alpha _s)^2}{w(a_A)w(a_B)}}G_{a_A/A}(x_A,\mu _f^2)G_{a_B/B}(x_B,\mu _f^2)D_{h/a_1}(z_1,M_f^2)`$
$`\times `$ $`{\displaystyle \frac{\alpha _s}{2\pi }}\left[A_0^v+A_0^s+A_0^{coll}\right]dx_Adx_Bdz_1d\mathrm{\Gamma }_2.`$
The three-body contribution, now evaluated in four dimensions, was given in Eq. (174) where now the soft and collinear regions of phase space are excluded.

As in the previous examples, the structure of the final result is two finite contributions, both of which depend on the soft and collinear cutoffs – one explicitly and one through the boundaries imposed on the three-body phase space. However, when both contributions are added while calculating an observable quantity, all dependence on the cutoffs cancels when sufficiantly small values of the cutoffs are used.

## IV Discussion and Conclusions

A technique for performing next-to-leading-logarithm calculations using Monte Carlo techniques was described in detail. The method uses two cutoff parameters which serve to separate the regions of phase space containing the soft and collinear singularities from the non-singular regions. The main derivations for experimentally degenerate, tagged, and heavy quark final states were given, as was a discussion of initial state factorization. We provided five illustrative examples applying the method.

The first example was that of the QCD corrections to electron-positron annihilation into a massive quark pair. The quark mass serves to regulate any would-be final state collinear singularities. The final state soft singular region is delineated using one cutoff. The second example was that of QCD corrections to electron-positron annihilation into a massless quark pair. In this case final state soft and collinear singularities are encountered. Both soft and collinear cutoffs are therefore required. They should be chosen such that $`\delta _c\delta _s`$. The third example of inclusive photon production in hadronic final states of electron-positron annihilation was presented, illustrating the use of fragmentation functions. Finally, the QCD corrections to lepton pair and single particle production in hadron-hadron collisions were given. These examples include both initial and final state soft and collinear singularities. The use of scale dependent parton distribution and fragmentation functions was explained.

The Monte Carlo results, integrated to give an inclusive cross section, were shown to be in complete agreement with those available in the literature. This is not the end of the utility of the method, but only the beginning. Given the full access to the parton four vectors and corresponding weights, we are free to combine them in any way that is consistent with an infrared-safe measurement function, which may include a jet finding algorithm and experimental cuts.

The method has been applied to a wide range of hard scattering processes and it has been found to be both simple to implement and numerically robust.

###### Acknowledgements.

We thank our collaborators Lewis Bergmann, Howie Baer, Jim Ohnemus, Bob Bailey, and Uli Baur, who have, over the years, applied and shaped the method presented herein, and Jack Smith for comments on the text. This work was supported in part by the U.S. Department of Energy, High Energy Physics Division, under contracts DE-FG02-97ER41022 and W-31-109-Eng-38.

## Appendix A Comparison with other methods

The essential difference between the phase space slicing and the subtraction methods may be gleaned from the following simple example ks . Consider an integral to be calculated:

$$I=\underset{ϵ0^+}{lim}\left\{_0^1\frac{dx}{x}x^ϵF(x)\frac{1}{ϵ}F(0)\right\},$$
(188)
where $`F(x)`$ is a known but complicated function related to a two-to-three body matrix element. The variable $`x`$ represents either the energy of an emitted gluon or the angle between two massless partons. In a traditional fully or single particle inclusive calculation the integral $`I`$ would be performed completely analytically.

In the subtraction method one simply adds and subtracts $`F(0)`$ under the integral sign.

$`I`$ $`=`$ $`\underset{ϵ0^+}{lim}\left\{{\displaystyle _0^1}{\displaystyle \frac{dx}{x}}x^ϵ\left[F(x)F(0)+F(0)\right]{\displaystyle \frac{1}{ϵ}}F(0)\right\}`$ (189)
$`=`$ $`{\displaystyle _0^1}{\displaystyle \frac{dx}{x}}\left[F(x)F(0)\right],`$
giving a finite and numerically calculable result. No approximations are made, however in any numerical implementation there will necessarily be a lower limit related to machine precision below which the integral must be cutoff. This is not a problem in practice.

In the phase space slicing method, the integration region is divided into two parts $`0<x<\delta `$ and $`\delta <x<1`$ with $`\delta 1`$. A Maclaurin expansion of $`F(x)`$ yields

$`I`$ $`=`$ $`\underset{ϵ0^+}{lim}\left\{{\displaystyle _0^\delta }{\displaystyle \frac{dx}{x}}x^ϵF(x)+{\displaystyle _\delta ^1}{\displaystyle \frac{dx}{x}}x^ϵF(x){\displaystyle \frac{1}{ϵ}}F(0)\right\}`$ (190)
$`=`$ $`{\displaystyle _\delta ^1}{\displaystyle \frac{dx}{x}}F(x)+F(0)\mathrm{ln}\delta +𝒪(\delta ).`$
Clearly, the parameter $`\delta `$ must be chosen small enough so that the term linear in $`\delta `$ may be neglected. At the same time it must not be so small as to spoil the numerical convergence of the first term.

## Appendix B Soft Integrals

In evaluating the soft integrals we encounter angular integrals which may be written in the form

$$I_n^{(k,l)}=_0^\pi 𝑑\theta _1\mathrm{sin}^{n3}\theta _1_0^\pi 𝑑\theta _2\mathrm{sin}^{n4}\theta _2\frac{(a+b\mathrm{cos}\theta _1)^k}{(A+B\mathrm{cos}\theta _1+C\mathrm{sin}\theta _1\mathrm{cos}\theta _2)^l}.$$
(191)
A large collection of these appear in the appendix of been . Others may be found in the appendix of hs  or else computed as explained in willy . Here we collect together the results covering most of the cases encountered using the two cutoff slicing method. The first two are from hs  with $`A^2B^2+C^2`$

$`I_n^{(0,1)}`$ $`=`$ $`{\displaystyle \frac{\pi }{\sqrt{B^2+C^2}}}\{\mathrm{ln}\left({\displaystyle \frac{A+\sqrt{B^2+C^2}}{A\sqrt{B^2+C^2}}}\right)`$
$`(n4)[\mathrm{Li}_2\left({\displaystyle \frac{2\sqrt{B^2+C^2}}{A+\sqrt{B^2+C^2}}}\right)+{\displaystyle \frac{1}{4}}\mathrm{ln}^2\left({\displaystyle \frac{A+\sqrt{B^2+C^2}}{A\sqrt{B^2+C^2}}}\right)]\}`$
$$I_n^{(0,2)}=\frac{2\pi }{A^2B^2C^2}\left[1\frac{1}{2}(n4)\frac{A}{\sqrt{B^2+C^2}}\mathrm{ln}\left(\frac{A+\sqrt{B^2+C^2}}{A\sqrt{B^2+C^2}}\right)\right],$$
(193)
where we drop $`O((n4)^2)`$ terms. The second two are from been  with $`b=a`$. If $`A^2=B^2+C^2`$

$$I_n^{(1,1)}=2\pi \frac{1}{aA}\frac{1}{n4}\left(\frac{A+B}{2A}\right)^{n/23}\left[1+\frac{1}{4}(n4)^2\mathrm{Li}_2\left(\frac{AB}{2A}\right)\right],$$
(194)
whereas if $`A^2B^2+C^2`$

$`I_n^{(1,1)}`$ $`=`$ $`{\displaystyle \frac{\pi }{a(A+B)}}\{{\displaystyle \frac{2}{n4}}+\mathrm{ln}\left[{\displaystyle \frac{(A+B)^2}{A^2B^2C^2}}\right]`$ (195)
$`+`$ $`{\displaystyle \frac{1}{2}}(n4)[\mathrm{ln}^2\left({\displaystyle \frac{A\sqrt{B^2+C^2}}{A+B}}\right){\displaystyle \frac{1}{2}}\mathrm{ln}^2\left({\displaystyle \frac{A+\sqrt{B^2+C^2}}{A\sqrt{B^2+C^2}}}\right)`$
$`+`$ $`2\mathrm{Li}_2({\displaystyle \frac{B+\sqrt{B^2+C^2}}{A\sqrt{B^2+C^2}}})2\mathrm{Li}_2\left({\displaystyle \frac{B\sqrt{B^2+C^2}}{A+B}}\right)]\},`$
again dropping $`O((n4)^2)`$ terms in the second of these. The dilogarithm function $`\mathrm{Li}_2(x)`$ is defined in dilog  and numerous useful properties are summarized in dd .

## Appendix C Recovering the $`𝒪(\delta _c/\delta _s)`$ terms

In this appendix we integrate the $`P_{qq}(z,ϵ)`$ splitting kernel over the hard-collinear portion of phase space for the case of a 45 singularity as it pertains to the discussion given at the end of Sec. III.2. Recall that this region is defined by

$`\mathrm{hard}`$ $`:`$ $`\delta _s{\displaystyle \frac{\sqrt{s_{12}}}{2}}E_5{\displaystyle \frac{\sqrt{s_{12}}}{2}}`$
$`\mathrm{collinear}`$ $`:`$ $`0s_{45}\delta _cs_{12}.`$ (196)
From Eq. (29) we have $`s_{34}=(p_3+p_4)^2=2p_3p_4(2p_3p_{45})z`$ and $`s_{12}=(p_3+p_{45})^2s_{45}+2p_3p_{45}`$ which together yield $`s_{34}z(s_{12}s_{45})`$. Using $`E_5=(s_{12}s_{34})/2\sqrt{s_{12}}`$ the hard condition becomes

$$0z\frac{1\delta _s}{1s_{45}/s_{12}}.$$
(197)
The approximation made in Sec. II.2.1 was to set $`s_{45}=0`$ in the denominator, in light of the collinear condition. This resulted in a decoupling of the $`z`$ and $`s_{45}`$ integration limits in Eq. (36). Relaxing the $`s_{45}=0`$ approximation gives rise to terms $`𝒪(\delta _c/\delta _s)`$ as now described.

Keeping the $`s_{45}`$ dependence, the required integral is

$$I=_0^{\delta _cs_{12}}\frac{ds_{45}}{s_{45}}\left(\frac{s_{45}}{s_{12}}\right)^ϵ_0^{\frac{1\delta _s}{1s_{45}/s_{12}}}𝑑z[z(1z)]^ϵP_{qq}(z,ϵ).$$
(198)
We may expand $`P_{qq}(z,ϵ)`$ about $`ϵ=0`$ and make a change of variables $`u=s_{45}/s_{12}`$ giving

$$I=C_F_0^{\delta _c}u^{1ϵ}F(u)𝑑u,$$
(199)
with

$$F(u)=_0^{\frac{1\delta _s}{1u}}𝑑z\left\{\frac{1+z^2}{1z}ϵ(1z)ϵ\mathrm{ln}\left[z(1z)\right]\frac{1+z^2}{1z}\right\}.$$
(200)
$`F(u)`$ may be evaluated with the help of

$`{\displaystyle _0^a}𝑑z{\displaystyle \frac{1+z^2}{1z}}`$ $`=`$ $`a\left(1+{\displaystyle \frac{a}{2}}\right)2\mathrm{ln}(1a)`$ (201)
$`{\displaystyle _0^a}𝑑z(1z)`$ $`=`$ $`a\left(1{\displaystyle \frac{a}{2}}\right)`$ (202)
$`{\displaystyle _0^a}𝑑z{\displaystyle \frac{\mathrm{ln}z}{1z}}`$ $`=`$ $`\mathrm{Li}_2(1a){\displaystyle \frac{\pi ^2}{6}}`$ (203)
$`{\displaystyle _0^a}𝑑z{\displaystyle \frac{\mathrm{ln}(1z)}{1z}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}^2(1a)`$ (204)
$`{\displaystyle _0^a}𝑑z{\displaystyle \frac{z^2\mathrm{ln}z}{1z}}`$ $`=`$ $`a\left(1+{\displaystyle \frac{a}{4}}\right)a\left(1+{\displaystyle \frac{a}{2}}\right)\mathrm{ln}a+\mathrm{Li}_2(1a){\displaystyle \frac{\pi ^2}{6}}`$ (205)
$`{\displaystyle _0^a}𝑑z{\displaystyle \frac{z^2\mathrm{ln}(1z)}{1z}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}a\left(1+{\displaystyle \frac{a}{6}}\right)+\mathrm{ln}(1a)\left[{\displaystyle \frac{3}{2}}a\left(1+{\displaystyle \frac{a}{2}}\right){\displaystyle \frac{1}{2}}\mathrm{ln}(1a)\right].`$ (206)
The resulting terms in $`F(u)`$ may be integrated over $`u`$ using

$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\left({\displaystyle \frac{1\delta _s}{1u}}\right)^i`$ $`=`$ $`\left({\displaystyle \frac{1}{ϵ}}+\mathrm{ln}\delta _c\right);i=0,1,2`$ (207)
$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\mathrm{ln}\left(1{\displaystyle \frac{1\delta _s}{1u}}\right)`$ $`=`$ $`\left({\displaystyle \frac{1}{ϵ}}+\mathrm{ln}\delta _c\right)\mathrm{ln}\delta _s\mathrm{Li}_2(\delta _c/\delta _s),`$ (208)
for the terms multiplied by $`𝒪(ϵ^0)`$ in Eq. (200) and

$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\mathrm{ln}^2\left(1{\displaystyle \frac{1\delta _s}{1u}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\mathrm{ln}^2\delta _s`$ (209)
$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\left({\displaystyle \frac{1\delta _s}{1u}}\right)^i\mathrm{ln}\left(1{\displaystyle \frac{1\delta _s}{1u}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\mathrm{ln}\delta _s;i=1,2`$ (210)
$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\left({\displaystyle \frac{1\delta _s}{1u}}\right)^i\mathrm{ln}\left({\displaystyle \frac{1\delta _s}{1u}}\right)`$ $`=`$ $`0;i=1,2`$ (211)
$`{\displaystyle _0^{\delta _c}}𝑑uu^{1ϵ}\mathrm{Li}_2\left(1{\displaystyle \frac{1\delta _s}{1u}}\right)`$ $`=`$ $`0,`$ (212)
for the terms multiplied by $`𝒪(ϵ^1)`$ in Eq. (200). Terms containing or leading to contributions of $`𝒪(\delta _c)`$ or $`𝒪(\delta _s)`$ have been dropped. Taking the coefficients of $`1/ϵ`$ and $`ϵ^0`$ gives the desired result

$`A_1^{qqg}`$ $`=`$ $`C_F\left(3/2+2\mathrm{ln}\delta _s\right)`$ (213)
$`A_0^{qqg}`$ $`=`$ $`C_F\left[7/2\pi ^2/3\mathrm{ln}^2\delta _s+2\mathrm{Li}_2(\delta _c/\delta _s)\mathrm{ln}\delta _c\left(3/2+2\mathrm{ln}\delta _s\right)\right].`$ (214)
The second equation is identically Eq. (47) with the addition of the advertised $`\mathrm{Li}_2(\delta _c/\delta _s)`$ term. A similar analysis may be performed for the $`P_{gg}`$ splitting case with the same result: $`\mathrm{ln}\delta _c\mathrm{ln}\delta _s\mathrm{ln}\delta _c\mathrm{ln}\delta _s\mathrm{Li}_2(\delta _c/\delta _s)`$.

## Appendix D Improving convergence of tilde terms

We want to demonstrate how the numerical convergence of the $`\stackrel{~}{D}`$ and $`\stackrel{~}{G}`$ functions may be improved. To this end consider the integral

$$F_i=_x^{1\delta _s}\frac{dy}{y}G_j(x/y,\mu )\stackrel{~}{P}_{ij}(y),$$
(215)
with

$$\stackrel{~}{P}_{ij}(y)=P_{ij}(y)\mathrm{ln}\left(\delta _c\frac{1y}{y}\frac{s_{12}}{\mu _f^2}\right)P_{ij}^{}(y).$$
(216)
Here a logarithm of $`\delta _s`$ is numerically being built up. Convergence will be improved if we rewrite the result in a form where the logarithmic dependence on $`\delta _s`$ is manifest. To do so, use the fact that

$$\underset{y1}{lim}\left[(1y)P_{ij}(y)\right]=2C_i\delta _{ij},$$
(217)
with $`C_g=N`$ and $`C_q=C_F`$. Now add and subtract the leading singular piece under the integral sign:

$`F_i`$ $`=`$ $`{\displaystyle _x^{1\delta _s}}{\displaystyle \frac{dy}{y}}[G_j(x/y,\mu )P_{ij}(y)\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)G_j(x/y,\mu )P_{ij}^{}(y)`$ (218)
$``$ $`G_j(x,\mu ){\displaystyle \frac{2C_i\delta _{ij}}{1y}}\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)+G_j(x,\mu ){\displaystyle \frac{2C_i\delta _{ij}}{1y}}\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)].`$
Regrouping terms gives

$`F_i`$ $`=`$ $`{\displaystyle _x^{1\delta _s}}{\displaystyle \frac{dy}{y}}\{[G_j(x/y,\mu )P_{ij}(y)G_j(x,\mu ){\displaystyle \frac{2C_i\delta _{ij}}{1y}}]\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)`$ (219)
$``$ $`G_j(x/y,\mu )P_{ij}^{}(y)+G_j(x,\mu ){\displaystyle \frac{2C_i\delta _{ij}}{1y}}\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)\}.`$
The last term may be evaluated with the help of

$$_x^{1\delta _s}\frac{dy}{y(1y)}\mathrm{ln}\left(a\frac{1y}{y}\right)=\frac{1}{2}\mathrm{ln}^2\left(\frac{a\delta _s}{1\delta _s}\right)+\frac{1}{2}\mathrm{ln}^2\left[\frac{a(1x)}{x}\right].$$
(220)
The final desired expression is

$`F_i`$ $`=`$ $`C_i\delta _{ij}G_j(x,\mu )\left[\mathrm{ln}^2\left(\delta _c{\displaystyle \frac{s_{12}}{\mu _f^2}}{\displaystyle \frac{1x}{x}}\right)\mathrm{ln}^2\left(\delta _c{\displaystyle \frac{s_{12}}{\mu _f^2}}{\displaystyle \frac{\delta _s}{1\delta _s}}\right)\right]`$ (221)
$`+`$ $`{\displaystyle _x^{1\delta _s}}{\displaystyle \frac{dy}{y}}\{[G_j(x/y,\mu )P_{ij}(y)G_j(x,\mu ){\displaystyle \frac{2C_i\delta _{ij}}{1y}}]\mathrm{ln}\left(\delta _c{\displaystyle \frac{1y}{y}}{\displaystyle \frac{s_{12}}{\mu _f^2}}\right)`$
$`G_j(x/y,\mu )P_{ij}^{}(y)\}.`$
The $`\mathrm{ln}\delta _s`$ is now evident in the first term, and absent from the second integral term. Numerical convergence will therefore be greatly improved.