Summing threshold logs in a parton shower

When parton distributions are falling steeply as the momentum fractions of the partons increases, there are effects that occur at each order in $\alpha_s$ that combine to affect hard scattering cross sections and need to be summed. We show how to accomplish this in a leading approximation in the context of a parton shower Monte Carlo event generator.


Introduction
One often wants to calculate the cross section for a hard process in hadron-hadron collisions under the circumstance that one or both of the required momentum fractions η a and η b for the initial state partons is large. Then the corresponding parton distribution functions f a/A (η a , µ 2 ) or f b/B (η b , µ 2 ) will be steeply falling as a function η a and η b . In this case, the cross section calculated beyond the leading order is enhanced compared to the leading order cross section. The effect is said to be due to "threshold logarithms." 1 The amount of enhancement is proportional to how fast f a/A (η a , µ 2 ) or f b/B (η b , µ 2 ) is falling. The threshold enhancement grows as η a and η b increase, but η a and η b do not have to be close to 1 for the threshold logarithms to start to be significant. One can count η a > 0.1 and η b > 0.1 as being part of the threshold region.
Even before the study of threshold logarithms began, Sjöstrand [48] and Marchesini and Webber [49] had introduced parton shower event generators that were important in the mid 1980s and whose direct descendants are still essential tools today [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65]. These programs are based on iterating parton splittings derived from perturbative quantum chromodynamics (QCD). Because of this iterative structure based on the soft and collinear singularities of the theory, they sum many sorts of logarithms in a natural way. For instance, parton shower event generators can sum logarithms of k 2 T /M 2 Z , where k T is the transverse momentum of a Z-boson produced in the Drell-Yan process [66]. Also, the analysis of ref. [4] connects the summation of threshold logarithms to the structure of a parton shower.
Current standard parton shower event generators do not sum threshold logarithms. However, it is possible to rearrange them so that they do sum the leading threshold logarithms. In this paper, we show how to do that.
Our analysis is based on our own parton shower event generator, Deductor [67][68][69][70][71][72][73][74]. We make this definite choice because we need to be specific with respect to shower kinematics, the choice of shower ordering variable, splitting functions, and the treatment of color. For this paper, we use a new version, 2.0.2, of Deductor [75].
Even though our analysis uses definitions specific to Deductor, we think that it will be apparent that the results could be adapted to other hardness based parton shower generators, such as Pythia [65] and Sherpa [61]. With somewhat more restructuring, it should be possible to adapt the results to the angle ordered parton shower generator Herwig [60].
We present some of the features of the theoretical framework in section 2. Then, before we begin the analysis, we exhibit in section 3 some of the main features of the result and of where this result comes from. Inevitably, we leave a lot out. The details are presented later in the paper. Section 4 discusses the role of parton distribution functions. Section 5 analyzes the Sudakov exponent conventionally used in parton shower algorithms. Section 6 presents a revised Sudakov exponent to be used if one wants to include threshold logarithms. Sections 7 and 8 analyze the exponential that contains the threshold logarithms. Section 9 presents a comparison of the leading behavior of the threshold factor in this paper to the factor in previous treatments of threshold logarithms. Section 10 presents numerical results obtained with the Deductor parton shower. There are three concluding sections. In section 11, we summarize the main steps of the analysis that has been presented, leaving out the supporting arguments. The reader may wish to read section 11 first, then return to it after absorbing the supporting arguments. In section 12, we describe some of the choices available in running Deductor. In section 13, we say something about the advantages and disadvantages of putting threshold enhancements into a parton shower. We have tried to keep the main presentation relatively short by relegating calculational details to three appendices.

The Deductor framework
Our analysis is based on the parton shower event generator, Deductor [67][68][69][70][71][72][73][74], which uses specific choices with respect to shower kinematics, the shower ordering variable, the parton splitting functions, and the treatment of color. In this section, we outline some of these choices that play a role in the analysis of this paper.
An exact treatment of leading threshold logarithms requires an exact treatment of color, which is available in the general formalism of Refs. [67][68][69][70][71][72][73][74]. The exact color treatment is not implemented in the code of Deductor. Rather, we are able to use only an approximation, the leading-color-plus (LC+) approximation. However, in this paper we develop the theory using exact color. Then the LC+ approximation consists of simply dropping some terms.
In Deductor, we order splittings according to decreasing values of a hardness parameter Λ 2 . The hardness parameter is based on virtuality. For massless partons, the definition is (2.1) Here the mother parton in a final state splitting has momentum p l and the daughters have momentap l andp m+1 . For an initial state splitting in hadron A, the mother parton has momentum p a , the new initial state parton has momentump a and the final state parton created in the splitting has momentump m+1 . 2 For hadron B, we replace "a" → "b." We denote by Q 0 a fixed vector equal to the total momentum of all of the final state partons just after the hard scattering that initiates the shower. We often use a dimensionless "shower time" variable given by t = log(Q 2 0 /Λ 2 ) .

(2.2)
Then t increases as the shower develops. One could use some other hardness parameter to order the shower. For instance, various measures of the transverse momentum in a splitting are popular choices. In this paper, we use Λ 2 because that is what we have implemented in Deductor.
In an initial state splitting, parton distribution functions enter the splitting functions. Then we need a scale parameter for the parton distribution functions. We use the virtuality of the splitting: (2.3) 2 Here pa is the momentum of the mother parton in the sense of shower evolution that moves from hard interactions to softer interactions. With initial state splittings, the shower development moves backwards in physical time.
Similarly, for a final state splitting of parton l, the virtuality parameter is µ 2 l (t) = 2p l · Q 0 e −t . (2.4) It proves convenient to use a dimensionless virtuality variable y. For an initial state splitting, where Q = p a + p b is the total momentum of the final state particles in the shower state before the splitting. We collect these and other notations used throughout the paper in appendix A.
Deductor 1.0 uses non-zero charm and bottom quark masses, both in the final and the initial states. However, we have simplified the code in Deductor 2.0 [75] that we use here by setting the charm and bottom quark masses to zero in the parton kinematics (although charm and bottom quark thresholds in the evolution of parton distribution functions remain). This is evidently an approximation, but this approximation does not matter for the physics investigated in this paper, the threshold logarithms at very large parton collision energies.
The general formalism includes parton spins, but the current implementation in Deductor simply eliminates spin by averaging over the mother parton spin before each splitting and summing over the daughter parton spins. In this paper we keep the analysis as simple as we can by doing the same spin averaging, so that spin dynamics does not play a role.
3 Some main features of the threshold enhancement As in our previous work [67][68][69][70][71][72][73][74], we find it useful to think of parton shower evolution as describing the evolution of a statistical state vector ρ(t) . The shower time t represents the hardness scale of the current state. As t increases, parton emissions get softer and softer. At shower time t, ρ(t) represents the distribution of partonic states in an ensemble of runs of the Monte Carlo style program. Without color or spin, ρ(t) would simply represent a probability distribution in the number m + 2 of partons (m final state partons and two initial state partons), their flavors f , and their momenta p. With the inclusion of color, ρ(t) is a little more complicated: it represents the probability distribution of parton flavors and momenta and the density matrix in the quantum color space. As stated in the introduction, we ignore spin in this paper. If spin were included, ρ(t) would represent the density matrix in the combined color and spin space.
In the notation we use here, a measurement of some property F of the partonic system, like the number of jets with transverse momenta bigger than some given value, is represented by a vector F . The average cross section corresponding to measurement F at time t is F ρ(t) . Of particular importance is the completely inclusive measurement, which we denote by 1 . Thus the total cross section represented by state ρ(t) is 1 ρ(t) . To compute 1 ρ(t) , we sum over the number of partons and their flavors, integrate over the parton momenta, and take the trace of the density matrix in color space.
We represent the evolution of the shower by the equation Here H I (t) represents real parton emissions and increases the number of final state partons by one. It contains splitting functions derived from QCD and a ratio of parton distribution functions. The operator S(t) leaves the number of partons, their momenta, and their flavors unchanged. It is, however, an operator on the color space. The operator S(t) contains terms representing two effects. It contains a term derived from one loop virtual graphs. It also contains a term proportional to the derivative with respect to scale of the parton distribution functions. This term is needed because parton distribution functions appear at the hard interaction and then (with a different scale) at each parton splitting. Thus one needs a term in S(t) that cancels the parton distributions at the old, harder, scale and inserts parton distributions at the new, softer, scale. It is useful to solve the shower evolution equation (3.1) in the form Here N S (τ, t 0 ) is the no-splitting operator, given by a time-ordered exponential of the integral of S(τ ),  That is, we define V(t) from H I (t) so that 1 V(t) = 1 H I (t) . (3.9) With this definition, the color trace of the revised Sudakov factor, represents the probability for the parton system not to split between shower time t 1 and shower time t 2 . It is standard to use Sudakov factors that are exponentials of V(t) defined by eq. (3.7) to define a parton shower. That is the choice in Deductor 1.0. When one constructs a parton shower using V(t), the cross section associated with ρ(t) is what it was at the hard interaction, say t = 0. That is to say, the total cross section is the Born cross section. What the parton shower does is to distribute the fixed cross section into cross sections for the different final states that the starting partons could evolve into.
The idea of the summation of threshold logarithms is that the total cross section is not just the Born cross section. Rather it contains corrections from higher orders of perturbation theory. To see this effect, we should use S(t) instead of V(t). When we do that, we have a Sudakov factor that we can write as (3.11) Compared to the standard Sudakov factor made using V(τ ), there is then an extra term in the exponent, [V(τ ) − S(τ )]. The most important parts of this are those associated with initial state evolution, There is a good deal of calculation needed to find [V a (τ ) − S a (τ )]. However, the most important term in the result is pretty simple to understand. For a partonic state with m + 2 partons with momenta p, flavors f , and colors c , c, we find that the contribution is (3.12) The omitted terms here need not concern us at the moment. The function P aa (z) is the flavor conserving part of the evolution kernel for the parton distributions f a/A (η a , µ 2 ). The factors [1 ⊗ 1] represent unit operators on the color space. The index a represents the flavor of parton "a" in the state {p, f, c , c} m . The factors C a are color factors, C F for quark or antiquark flavors, C A for gluons. The parameter y, defined in eq. (2.5), is a dimensionless virtuality variable for a splitting at shower time t. For large t, the splitting has small virtuality and y 1. The integration variable z is the momentum fraction in the splitting: parton "a" carries momentum fraction η a before the splitting and the new parton "a" carries momentum fraction η a /z after the splitting.
The first term in eq. (3.12) comes from V a (t), which is given in eq. (5.6) and involves parton distribution functions and the splitting functions that describe real emissions. The second term comes from S a (t), which is given in eq. (6.18) and involves the evolution of the parton distribution functions together with contributions from virtual graphs. The full V a (t) − S a (t) is given in eq. (7.6).
In the first term in eq. (3.12), the upper limit z < 1/(1 + y) comes from the kinematics of parton splitting: for a finite virtuality, (1 − z) cannot be too small. See eq. (A.7). In the second term, the upper limit z < 1 comes from the upper limit in the evolution equation for parton distribution functions. We see that the two terms are almost identical. They differ only in the upper limits of the z-integrations. Adding the two terms, we have Only a small integration range, corresponding to soft gluon emission, remains. In this range, we can use This gives We integrate over a small range of z near z = 1. There is a 1/(1 − z) singularity, the "threshold singularity." The ratio of parton distribution functions is 1 at z = 1. Thus the integral is not singular. But the integral is large if f a/A (η a , µ 2 ) is a fast varying function of η a . When η a is not small, say η a > 0.1, f a/A (η a , µ 2 ) is in fact a fast varying function of η a . For that reason, we retain this contribution. The difference [V a (t) − S a (t)] appears in the exponent of the Sudakov factor, so by keeping this contribution, we sum the effects that come from the threshold singularity. This "shower" version of the summation of threshold logs looks rather different from the direct QCD or the SCET descriptions of the same physics. We will discuss the connections in section 9.
With this introduction to provide hints about where we are going, we now proceed to a detailed examination of how the leading threshold logarithms can be summed within a parton shower event generator, at least after we apply the LC+ approximation [70] to simplify the color structure.

The parton distribution functions
As was noted already in ref. [1], the precise definition of parton distributions matters for the question of threshold logarithms. In this paper, we are dealing with a parton shower algorithm. As we will see in more detail in the following sections, it is important that the behavior of the parton distribution functions match the structure of the parton shower algorithm. We discuss the needed definitions in this section.

Parton evolution equations
Begin with the MS parton distribution functions with one loop evolution. Call the scale variable µ 2 ⊥ . We have (4.1) In the case of a bottom or charm quark, evolution of f a/A (η a , µ 2 ⊥ ) occurs only for µ 2 ⊥ > m 2 ⊥ (a), where m ⊥ (a) is the quark mass. Below this scale, f a/A (η a , µ 2 ⊥ ) = 0. 3 For a light quark or gluon, we define m 2 ⊥ (f ) = m 2 ⊥ (start) to be a starting scale for shower evolution and for parton distribution function evolution. In Deductor, we take m 2 ⊥ (start) near 1 GeV 2 . In eq. (4.1), P aâ (z) are the familiar DGLAP kernels and C a and γ a are the standard flavor dependent constants recorded in eqs. (A.9) and (A.10). We take the argument of α s to be That is, we use the MS scale µ 2 ⊥ for α s in the evolution equation and we insert an extra factor λ R . (In λ R , the number n f of active flavors depends on µ 2 ⊥ .) The factor λ R in the argument of α s is, strictly speaking, beyond the order of perturbation theory that we control in a leading order shower, but it is helpful in generating certain next-to-leading logarithms [51,66].
To make the structure of the z → 1 singularities in this equation clearer, we write where P reg aâ (z) is nonsingular as z → 1. Then we can write .

(4.4)
Now we look at the parton evolution that matches shower evolution in Deductor, which uses the virtuality based ordering variable Λ 2 defined in eq. (2.1). It thus needs parton distribution functions f a/A (η, µ 2 Λ ), where µ 2 Λ represents the virtuality in an initial state splitting, as in eq. (2.3). In the conventional MS parton distribution functions [76], the scale variable is the renormalization scale for the functions with MS renormalization. In MS renormalization at one loop, we have a dimensionally regulated integration, over the momentum fraction z and the transverse momentum k ⊥ of the splitting. The integration is ultraviolet divergent at = 0. The MS prescription is to subtract the 1/ pole along with certain conventional finite terms. This gives us a finite result in which the scale of the squared transverse momentum k 2 ⊥ is set by µ 2 ⊥ . This is not quite the same as setting the scale of the virtuality |k 2 | in the splitting. The two variables are related, 4 for small angle splittings, by k 2 , evolution for partons of flavor a stops. We want parton evolution and shower evolution to match. Thus we define With this definition of the scales, under under an infinitesimal change of µ 2 Λ we have . (4.7) Thus the one loop evolution equation for virtuality based parton distribution functions is almost the same as the evolution equation for MS parton distribution functions. We need theta functions that turn off the evolution for flavor a unless (1 − z)µ 2 Λ > m 2 ⊥ (a) and we should use (1 − z)λ R µ 2 Λ in the argument of α s . We simplify this a bit and use . (4.8) In the first term, there is no change from the MS evolution equation. In the remaining two terms, there is an integration over z and we identify µ 2 ⊥ with (1 − z)µ 2 Λ in the theta function that provides an infrared cutoff on µ 2 ⊥ . In the second term, which has a 1/(1 − z) singularity from soft gluon emission, we use (1 − z)λ R µ 2 Λ as the argument of α s , while in the third term, with no 1/(1 − z) singularity, we simplify the evolution by omitting the factor (1 − z) in the argument of α s .
With these choices, parton evolution matches the conventions that we use in the shower splitting kernels in Deductor. The parton distributions f a/A (η a , µ 2 Λ ) that we use are obtained by solving the evolution equation (4.8) using an MS parton distribution set as the initial condition at the starting scale m 2 ⊥ (start) for the shower. For the starting distributions, we use the CT14 NLO parton set [77].

Approximate analytic result
In section 9, we will present an analytical comparison of the results of this paper for the Drell-Yan cross section to standard analytical results in the leading log approximation. For this purpose, we will need an approximate analytical relation between the parton distribution functions f a/A (η a , µ 2 Λ ) obtained by solving eq. (4.8) and the MS parton distribution functions f MS a/A (η a , µ 2 ⊥ ) obtained by solving eq. (4.4). Our aim is to understand just the leading contributions to threshold behavior, so we do not specify the argument of α s and we make some approximations that correspond to including only the effect of soft gluon emissions.
We write eqs. (4.8) and (4.4) as equations for the Mellin moments of the parton distributions,f We introduce a parameter λ with λ = 0 corresponding to MS evolution, eq. (4.4), and λ = 1 corresponding to Λ 2 evolution, eq. (4.8). We denote by µ 2 0 the appropriate scale parameter, µ 2 ⊥ for λ = 0 or µ 2 Λ for λ = 1. Keeping only the leading singular terms in the kernel, we have (4.10) The solution of this is At the starting scale, m 2 ⊥ (start), we use the MS parton distributions (λ = 0) as a boundary value.
Now we can regard λ as a continuous variable. Then eq. (4.11) gives us a differential equation for the λ dependence off a/A (N, µ 2 0 , λ): (4.12) When we take the inverse Mellin transform of this, we find .

(4.17)
When does Z a (η a , µ 2 0 ) differ significantly from 1? There is a factor α s in the exponent, which generally makes the exponent small. However, the factor multiplying α s can be large when the parton ratio R is a fast varying function of z. We will use this result in section 9.4.

The probability preserving integrand V(t)
As outlined in section 3, we seek to calculate (with suitable approximations) the difference (V(τ ) − S(τ )) between the integrand of the Sudakov exponent V(τ ) that preserves probabilities as the shower evolves and the Sudakov integrand S(τ ) that is based on virtual diagrams and the evolution of the parton distribution functions.
We begin in this section with V(τ ) as defined in Deductor 1.0, but with all quark masses set to zero. We sketch the needed calculations in appendix B. We will use the results in section 7.
The operator V(τ ) is used to define the Sudakov factor, The exponent , when we take its trace in color, is the total probability to have a real parton splitting between shower times t 1 and t 2 . Thus N V (t 2 , t 1 ) is the probability not to have had a splitting. The operator V(τ ) has, in general, a non-trivial color structure. For that reason, the exponential function is ordered in shower time τ , as indicated by the T instruction. In Deductor, we apply an approximation, the LC+ approximation, that makes V(τ ) diagonal in color. Then the T ordering instruction is not needed. In this paper, we keep the full color structure for all operators, but then at the end we can apply the LC+ approximation. For many of our formulas for V(τ ), it is convenient to define energies and angles in a reference frame in which the total momentum Q of all of the final state partons has the form Q = (E Q , 0). Then p l is the momentum of parton l. The angle θ kl between two parton momenta is defined in this same reference frame. Thus It is convenient to define See appendix A, where C f l and γ f l are also defined.
In our analysis, we use a dimensionless virtuality variable y given by eqs. (2.5) and (A.3). The variable y is fixed by the shower time τ according to eq. (A.4). The structure of V(τ ) is rather complicated, but we simplify it by taking the leading behavior as y → 0. That is, we use the leading behavior of the splitting functions in the limits of soft or collinear splittings.
The parton shower in Deductor has an infrared cutoff: no splittings with a transverse momentum smaller than m ⊥ (start) ∼ 1 GeV are generated. In terms of y and the momentum fraction z in the splitting, the cutoff is y(1 − z) > m 2 ⊥ (start)/(2p a · Q) for an initial state splitting and yz(1 − z) > m 2 ⊥ (start)/(2p l · Q) for an final state splitting. However, in calculating the leading y 1 behavior of V(τ ), we assume that y is not so small that we need to be concerned with this infrared cutoff. In fact, we will find that very tiny values of y are not relevant for the threshold effects that we investigate. Thus we simply calculate V(τ ) with no infrared cutoff.
Some of the terms in V(τ ) depend on the angle θ kl between the splitting parton l and a parton k that forms part of a color dipole with parton l. This angle dependence arises because soft gluon emission from a color dipole depends on the angles among the partons. The angle θ kl is small when partons k and l are the daughter partons of a previous splitting that was nearly collinear. When this previous splitting was a final state splitting, ordering in Λ 2 for the new splitting of parton l implies a l y < 1 − cos θ kl . When a previous splitting that produced partons k and l was an initial state splitting with a small momentum fraction z kl , one can also have a l y > 1 − cos θ kl . This happens in the case of multi-regge kinematics, as discussed in section 5.4 of ref. [72]. This is the opposite kinematic regime from that of threshold logarithms, so we ignore this possibility in this paper. Thus we assume a l y 1 − cos θ kl in evaluating V(τ ). The terms in V(τ ) contain operators like [(T l · T k ) ⊗ 1] that act on the color space. This means that, in the ket part of the color density matrix, color generator T a acts on the color of parton l, color generator T a acts on the color of parton k, and we sum over the octet color index a. In this example, a unit operator acts on the bra part of the color density matrix. In manipulating the color operators, we use the identity k T k = 0 and the identity (T l · T l ) = C f l , where C q = C F for a quark or antiquark flavor q and C g = C A for a gluon.
At any stage in the shower, there are initial state partons "a" and "b" as well as final state partons l with l ∈ {1, . . . , m}. There is a term in V(τ ) for each parton: The limiting form for V l (τ ) for a final state parton, from eq. (B.35), is quite simple, × {p, f, c , c} m .
The limiting form for V a (τ ) for initial state parton "a" is not quite so simple because it involves parton distribution functions. From eq. (B.68), we have × {p, f, c , c} m . (5.6) The first two terms involve ratios of parton distribution functions together with subtractions. The ratios of parton distribution functions arise naturally in the Sudakov factors for initial state parton shower evolution from hard emissions to soft emissions [48]. Notice that there is an upper limit for the momentum fraction variable in the splitting: z < 1/(1 + y).
With the definition of z used in Deductor, z → 1 implies that the emitted parton has zero momentum, so that the virtuality variable y must vanish. Thus at finite y there must be a limit on z. The precise limit depends on the splitting kinematics used in Deductor. See eq. (A.7). In the second term, there is a non-trivial color factor and a function defined in eq. (B.64), where ψ ak = (1 − cos θ ak )/ 8(1 + cos θ ak ), as in eq. (A.13). For (1 − z) y/ψ ak , ∆ ak is approximately 1/(1 − z). However, for (1 − z) y/ψ ak , ∆ ak tends to zero faster than 1/(1 − z). Thus, when we integrate over z, the important integration region is y/(1 + y) < (1 − z) y/ψ ak provided that ψ ak 1. If ψ ak 1, there is no important integration region.
In the final two terms, there is no dependence on parton distribution functions and the z-integration has been performed. One term depends on θ ak , while the other term, with a trivial color factor, depends only on y.

The Sudakov integrand S(t)
In this section, we study the integrand S(t) for the Sudakov exponent. As in our discussion of V(t), we calculate without an infrared cutoff associated with the end of the shower at k 2 ⊥ = m 2 ⊥ (start) ∼ 1 GeV 2 . Our first task is to identify the parts of S(t) that come from parton evolution and from virtual graphs.

Decomposition of S(t)
The operator S(t) affects the evolution of ρ(t) through eq. (3.1), The color density matrix ρ(t) is defined to include the proper factors of parton distributions so that the color trace of ρ is a differential cross section [67]. It is related to the density matrix without parton distribution functions by where {p, f, c , c} m .
In eq. (6.4), is the perturbative splitting function, with no factors representing parton distribution functions. That is, the definition of H I (t) based on eq. (8.26) of ref. [67] contains explicit factors of ratios of parton distributions; to define H pert I (t), we remove these factors. We then interpret −S pert (t) in eq. (6.4), as giving the part of the evolution of ρ pert (t) that does not involve parton splitting. Thus we will calculate −S pert (t) using suitable approximations to virtual one loop Feynman graphs. We do that in appendix C.
Using eqs. (6.1), (6.2), (6.4), and (6.5), we see that the operator S(t) is related to S pert (t) by Acting on a statistical basis state, with a corresponding expression for λ F b . Using the parton evolution equation, this is To be more precise, we should use the evolution equation (4.8) appropriate for Λ 2 as the evolution variable. For the moment, it suffices to use this simple form: we ignore the theta function that enforces an infrared cutoff µ 2 ⊥ > m 2 ⊥ (a) and we do not specify the argument of α s .

Approach to calculating S pert (t)
We have argued that −S pert (t) should be calculated from virtual one loop Feynman diagrams. There are some choices for how to do that.
First, we should choose a gauge. Ultimately, we expect that Feynman gauge provides the most powerful method, especially if one wants a method that can be extended to higher orders of perturbation theory. However, in Feynman gauge, factorization of softer interactions from harder interactions does not work graph by graph. Gluons that are nearly collinear to external legs of a diagram or are very soft but that carry unphysical longitudinal polarization can couple to the interior of ultraviolet dominated subgraphs. These attachments can be treated after a sum over graphs by the use of the gauge invariance of the theory. However, if we use a physical gauge, only transversely polarized gluons can propagate over long distances. Then we do not need to sum over graphs. For this reason, in appendix C, we use Coulomb gauge.
Second, virtual graphs do not come with a definition of the shower time t, so to define −S pert (t) from virtual graphs, we need to provide a relation to t. Ultimately, we expect that the most powerful method is to consider (6.10) Here we integrate over the shower time with a lower cutoff at t. This integral is infrared divergent and we understand it to be regulated with dimensional regularization. Then we can obtain S pert (t) by differentiating with respect to t. This formulation gives us the possibility of incorporating the virtual graphs into a general theory that includes the real graphs and a systematic treatment of factorization and renormalization.
In this paper, we choose a method that is computationally a little involved, but relates the definition of the shower time in virtual graphs quite directly to that in the shower splitting functions. We recognize that the splitting functions arise simply from cut Feynman diagrams for real parton emissions in an approximation in which one is close to the infrared singular soft or collinear limits. There is a three dimensional integral over the momentum of the emitted parton, with a delta function that fixes this momentum as a function of the shower time.
For virtual diagrams, one can imitate the structure of the real emission diagrams. We begin with an integration over four components of the loop momentum q. Working in a reference frame in which the total momentum Q of the final state partons has only a time component, we can perform the integration over the energy q 0 in the loop. This leaves us with an integration d q over the three-momentum in the loop. Then we have an integrand with a similar structure to that of the real emission diagrams, but with i terms in the denominators modified.
For both real and virtual diagrams, we obtain integrands corresponding to the exchange of a soft gluon between partons l and k that form a color dipole. We partition each of these contributions into two, one with a collinear singularity when the gluon momentum q is collinear to the momentum of parton l and one with a collinear singularity when q is collinear to the momentum of parton k.
In the splitting function corresponding to a real emission diagram, when a parton with momentum q is emitted from a parton with velocity v l , there is an infrared singularity when | q | − q · v l → 0. This happens in the soft limit | q | → 0 or in the collinear limit that q becomes collinear with v l . Note that (| q | − q · v l ) is proportional to the virtuality associated with the emission: The shower time t associated with this splitting is defined by Here Q 0 is the total momentum of the final state partons at the start of the shower, as in section 2.
In the virtual diagrams, we find that the infrared singularities are controlled by factors of the form 1/(| q | − q · v l ). We simply insert a delta function that identifies this denominator with the shower time according to eq. (6.12). This leaves an integration over two dimensions, which we can arrange to have the form dz dφ to match the integrations in splitting functions.
With these manipulations, there is a term in S pert (t) for each parton in the shower at shower time t: In the following two subsections, we state our results for these terms. In each case, the results apply in the limit that the dimensionless virtuality variable y is small: y 1 and y (1 − cos θ kl ).

Virtual contributions for final state partons
We sketch the calculation for S pert × {p, f, c , c} m . (6.14) The constants C f and γ f are defined in appendix A. In assembling this result, we have used the color identities T l · T l = C f l and k T l · T k = 0. The most notable feature of eq. (6.14) is that S pert l (t) has contributions proportional to iπ. There is one such contribution for each final state parton k that can form a color dipole with parton l, but there is no contribution for k = a or k = b. These contributions persist no matter how small e −t is.

Virtual contributions for initial state partons
We sketch the calculation for S pert l (t) for initial state partons in appendix C. After combining the definitions (C.2), (C.
In assembling this result, we have used the color identities T a · T a = C a and k T a · T k = 0 and the kinematic identity E Q = 2| p a |.
As in the final state S pert l (t), there is an iπ contribution, this time associated with the color dipole formed by parton "a" and parton "b."

The complete S(t)
The complete operator S(t) is a sum of pieces associated with the individual partons, For the final state partons, we use simply S pert l (t) from eq. (6.14). For initial state parton "a", we use S pert a (t) from eq. (6.15) and add the contribution from parton evolution according to eq. (6.6), (6.7) and (6.9). For final state partons, this gives For initial state partons, there is a small simplification because terms proportional to γ a cancel, giving We will use these results in the following section.

The cross section changing exponent
The operator inserted between parton splittings, will preserve the Born cross section throughout the shower. This statement is exact in color as long as we use the exact V(τ ) and the exact splitting operator H I (t). In the simplest application, we use the leading color approximation or the LC+ approximation for both H I (t) and V(τ ). Then the Born cross section is also preserved in the shower. In a leading order shower, it would be better to use where S(t) is obtained from virtual one loop graphs and parton evolution as explained above. Then we have 3) The first term creates a cross section preserving shower, while the second term sums effects that change the cross section. We thus need the cross section changing integrand (V(τ ) − S(τ )) and its integral over τ . We assemble this from our previous results.

Final state partons
For the contribution from final state partons, we subtract eq. (6.17) for S l (τ ) from eq. (5.5) for V l (τ ). Almost everything cancels and we are left with × {p, f, c , c} m .
That is, only the "iπ" terms remain. We have supplied the virtuality of the splitting as the argument of α s , using the definition (2.4).

Initial state partons
For the contribution from initial state parton "a," we subtract eq. (6.18) for S a (τ ) from eq. (5.6) for V a (τ ). Quite a lot cancels and we are left with × {p, f, c , c} m .
There is an "iπ" term that does not cancel. There is a term with a non-trivial color factor that is proportional to the function ∆ ak (z, y) defined in eq. (5.7). Finally, there is a term with a trivial color factor [1 ⊗ 1]. This term arises from using a contribution from V a (t) that is integrated over 0 < z < 1/(1 + y) and subtracting a term with the same integrand from S a (t) that is integrated over 0 < z < 1. In the difference, we integrate over 1/(1 + y) < z < 1. For small y, this is a small range near z = 1.
We can integrate this over a range of shower times, t 1 < τ < t 2 . This gives the probability changing Sudakov exponent associated with going from a splitting at shower time t 1 to a later splitting at shower time t 2 . We use y = µ 2 a (τ )/Q 2 , where Q 2 is the square of the total momentum of the final state particles at shower time t 1 . Additionally, in the [1⊗1] term, we decompose the splitting function P aâ (z) into P aâ (z) = δ aâ 2zC a /(1−z)+P reg aâ (z) as in eq. (4.3). We keep the nonsingular P reg aâ (z) term because it can give a large contribution for certain flavor choices. This gives In the remainder of this section, we discuss the physics of the four terms in eq. (7.7), revise the ∆ ak term to better reflect the physics, implement an infrared cutoff, and specify the argument of α s in each term.

Integration ranges
We first consider the effective ranges of µ 2 and z that appear in each term in eq. (7.7). All but the last term in eq. (7.7) involve the ratio y = µ 2 /Q 2 for a potential splitting that might occur in hadron A between the previous splitting at shower time t 1 and the next shower time t 2 at which a splitting occurs. Suppose thatỹ =μ 2 /Q 2 is the corresponding hardness variable for a previous splitting in hadron A, perhaps the one at t 1 if it was in hadron A or perhaps an earlier splitting. How is µ 2 /Q 2 related toμ 2 /Q 2 ? The relation between µ 2 /Q 2 and the shower time for initial state splittings is given in eq. (A.4), As the shower progresses with increasing t, the momentum p b can stay the same or increase. Thus µ 2 /Q 2 decreases: We now examine the terms in eq. (7.7), beginning with the simplest term, the one proportional to iπ. This term can have important effects, but it does not change the inclusive probability 1 ρ(t) associated with the statistical state. Its effect continues for all µ 2 from the beginning of the shower to the end of the shower. Consider next the second term, proportional to P reg aâ (z). The integrand here has no z → 1 singularity. The integration over z covers the range 0 < 1 − z < µ 2 /Q 2 /(1 + µ 2 /Q 2 ). Thus this term contributes only near the start of the shower, when µ 2 /Q 2 ∼ 1. It does not contribute when µ 2 /Q 2 1. Now consider the first term. Here there is a factor 2C a /(1 − z), so we are sensitive to small (1 − z). We can understand the integration region, at least qualitatively, by approximating 5 the parton distribution functions as f a/A (η a , µ 2 ) ∝ η −N a . We consider the case that N is large, which happens when η a is bigger than about 0.1. Here "N is large" means something like N ∼ 3 in realistic cases. With this approximation, we have This factor in eq. (7.10) is 1 for (1 − z) 1/N and tends to zero for ( We also approximate 1 Thus the integral in the first term is roughly Evidently, this is a rather crude approximation, but it is instructive. The integral would be a small perturbative correction that we could simply ignore except that N is large. This gives a contribution to the exponent proportional to log 2 N . However, the range of µ 2 /Q 2 does not extend down to infinitesimal values. As soon as µ 2 /Q 2 < 1/N , there is no more contribution. That is, the threshold factor associated with this term in eq. (7.7) comes from the first few steps in shower evolution. The remaining term in eq. (7.7) has more complicated structure. It contains a sum over final state partons k and factor ∆ ak , defined in eq. (5.7). This factor depends on the angle θ ak the initial state parton "a" and parton k. The angle θ ak appears in the combination The factor 2 here is rather arbitrary.) For the parton factor in this term, we can use the rough approximation (7.11). We can also approximate the upper limit of the z-integration as (1 − z) < µ 2 /Q 2 . This gives us the integral (7.14) We will consider the behavior of this integral both for ψ ak ∼ 1 and for ψ ak 1. Let us first consider the case that ψ ak ∼ 1. Then we integrate (1 − z) over a range that is large when N is large, giving a log N . On the other hand, µ 2 /Q 2 is integrated over a finite range. Thus we have only a single power of log N . We note also that µ 2 /Q 2 is never smaller than a number of order 1/N . Now consider the case that ψ ak 1, supposing that µ 2 (t 2 ) → 0. Then the lower bound 2ψ ak (1 − z) on µ 2 /Q 2 in eq. (7.14) is very small and one may wonder if this leads to a large integral. When ψ ak 1, we have Then parton k makes a very small angle with respect to the beam axis. This can happen with a large probability when, late in the shower, parton k was emitted from the initial state parton from hadron A. The most important case to consider is that the latest real emission before the shower time interval under consideration was the emission of parton k. Then the upper endpoint of the µ 2 integration, µ 2 < µ 2 (t 1 ) corresponds to the bound (7.9), in whichμ 2 andQ 2 refer to the splitting at which parton k was created. We havẽ HereẼ a =zE a is the momentum of the mother parton of the previous splitting, while the energy of the of the emitted parton is approximately E k ≈ (1 −z)E a . Thus Since we define energies and angles in the Q = 0 frame, we have E a = E b . Thus The splitting function for the previous emission contained a ratio f a/A (η/z,μ 2 )/f a/A (η,μ 2 ) of parton distributions. Applying the approximation (7.11) to this ratio, we have (1 −z) < 1/N . Thus also,z is close to 1. This gives In the integrand of eq. (7.14), we have (1 − z)θ 2 ak /4 < µ 2 /Q 2 and 1/N < (1 − z). This gives When we combine the upper bound (7.19) with the lower bound (7.20), we see that µ 2 /Q 2 can vary only over a small range around θ 2 ak /(4N ). Similarly, (1 − z) varies only over a small range around 1/N .

Revised ∆ ak term
We conclude that there are no log N factors associated with the integration in the ∆ ak term eq. (7.7) in the case that ψ ak 1. The integral does have a finite contribution proportional to α s with no log N factors. However, a first order parton shower is not adequate to calculate this contribution accurately: the definition of the parton shower splitting functions incorporates the strong ordering condition µ 2 /Q 2 μ 2 /Q 2 and that condition is violated here. Ordinarily, the inclusion of an inaccurately calculated small perturbative correction to the cross section would be of little consequence. However, this correction can occur many times as the shower progresses, leading to a large, inaccurately calculated, correction. Thus, we eliminate this contribution from the ∆ ak term at small ψ ak by multiplying this term by θ(ψ ak > ψ min ), where the default value of ψ min is 10 −2 .

Infrared cutoff
Recall now from section 4.1 that the shower algorithm has an infrared cutoff that vetoes initial state splittings unless the transverse momentum in the splitting is above a minimum: Here m ⊥ (a) is the shower cutoff scale of order 1 GeV or the heavy quark mass in the case of a charm or bottom quark. Eq. (7.7) was derived with no infrared cutoff, but we can insert the cutoff by inserting a factor θ((1 − z)µ 2 > m 2 ⊥ (a)) into the integrations over z. This cutoff has negligible effect as long as µ 2 m ⊥ (a)Q. Since the whole integral is negligible unless µ 2 > Q 2 /N , we see that inserting the infrared cutoff has negligible effect as long as m ⊥ (a) < Q/N . This is the case in situations of phenomenological interest, in which Q is of order 1 TeV, or at least 100 GeV for Tevatron studies, and m ⊥ (a) is at most 5 GeV. Thus we make eq. (7.7) consistent with the rest of the shower algorithm and with our treatment of parton evolution by inserting a factor θ((1 − z)µ 2 > m 2 ⊥ (a)) into the z-integrations in Eq. (7.7). 6

Running coupling
Now we need to specify the argument of the running coupling α s . The simplest choice would be α s (µ 2 ). Instead, in the two soft-sensitive terms (proportional to 1/(1 − z) and ∆ ak ) we use λ R µ 2 ⊥ = λ R (1 − z)µ 2 as the argument of α s . Here λ R [51] is the constant defined in eq. (4.2). The factors (1 − z) and λ R in the argument of α s are, strictly speaking, beyond the order of perturbation theory that we control in a leading order shower, but it is helpful in generating next-to-leading logarithms for at least some inclusive observables [51,66]. The use of α s (λ R (1 − z)µ 2 ) inside the z integrations would create an artificial problem if we integrated down to (1 − z) = 0. However, with the cutoff (1 − z)µ 2 > m 2 ⊥ (a) to keep us out of the nonperturbative region, we do not encounter this problem. For the remaining two terms in eq. (7.7), we use simply α s (λ R µ 2 ).

Result
With these substitutions, we have × {p, f, c , c} m .
Some discussion of this result may be useful. In the first term here, there is a singular factor 1/(1 − z). As discussed above, the singularity is cancelled because the ratio of parton distribution functions approaches 1 as (1 − z) → 0. This constant will be large if the derivative of f a/A (η a , µ 2 a (t)) with respect to η a is large. Then there is a "threshold enhancement" of the cross section. Since the integration range in z disappears when µ 2 /Q 2 → 0, the important contribution comes from the region in which µ 2 /Q 2 is not too small.
In the second term, there is no singular factor 1/(1 − z). In standard treatments, there is no "threshold log." However, the ratio of parton distributions can be large, for instance for a = g andâ = u. For this reason, we retain this contribution.
The term in eq. (7.21) proportional to ∆ ak appears when there is a final state parton that can form a color dipole with the initial state parton. This does not happen in the starting configuration of the Drell-Yan process, before any final state partons have been emitted. However, the ∆ ak term does appear in the starting parton configuration for jet production. The ∆ ak term also appears for any hard process once one or more final state partons have been emitted by initial state radiation. The color factor for this term is nontrivial. If we use the LC+ approximation instead of full color, the ∆ ak term contributes when parton k is color connected to parton "a." For the reasons given above in this section, we turn this term off when ψ ak is too small. Here ∆ ak (z, y) and ψ ak were defined in eqs. (B.64) and (B.59).
The last term in eq. (7.21) is proportional to iπ. This term is not associated with 1/(1 − z) singularities, but it is potentially important.

The cross section changing exponent in the LC+ approximation
The color structure of V(τ ) and of (V(τ ) − S(τ )) is non-trivial. However, the current version of Deductor uses the LC+ approximation. In this approximation, the operators are diagonal in color and thus commute with each other. Then Thus, within the LC+ approximation, we can generate the shower using V(τ ) and then, at each splitting, multiply by a numerical factor K(t 2 , t 1 ). The exponent in K is a sum over partons The contributions for final state partons l are given in eq. (7.5). The contribution for initial state parton "a" is given in eq. (7.21) and for parton "b" we simply have to substitute a ↔ b. With exact color, the operators in the exponent change the color state. However, in the LC+ approximation these operators are color diagonal. For a splitting of parton l with helper parton k, [T l · T k ⊗ 1] vanishes unless k is color connected to l in the ket state and, if k is color connected to l, equals C A /2 if parton l is a gluon and equals C F if parton l is a quark or antiquark. For [1 ⊗ T l · T k ], we have the same factors if k is color connected to l in the bra state. This rule applies for l and k being either initial state or final state parton indices.

Comparison to the standard summation
In this section, we consider the cross section dσ/(dQ 2 dY ) to produce a muon pair with squared momentum Q 2 and rapidity Y and compare our result to standard results. To do that, we need two manipulations, which are interesting in their own right.

The single power approximation for the Mellin transform
A parton distribution function can be expressed as an integral over its Mellin transform, Heref is the Mellin transform of f and is a function of the Mellin variable n = N + i ω.
The integration contour runs from N − i ∞ to N + i ∞, where N is chosen such that the contour runs to the right of any singularities off a/A (n, µ 2 ). In the standard method, it is not the parton distribution function that appears in an exponent, but rather the Mellin moment variable n. However, there is a simple method that allows us to compare the results of this paper to the standard results. We note that, with a reasonable model of the behavior of f a/A (η a , µ 2 ), its Mellin transformf a/A (n, µ 2 ) has a saddle point at some point n = N along the real axis. If we choose to let the integration contour in eq. (9.1) run through the saddle point, the integration will be dominated 7 by n ≈ N . (Cf. ref. [27] for a related use of the saddle point approximation.) Now, when addressing threshold summation, we encounter f a/A (η a /z, µ 2 ), where η a is the momentum fraction that appears in the Born cross section and we integrate over z.
Thus what enters our calculation is The integration over z is dominated by z near 1. Using the saddle point approximation, with the location N of the saddle point determined for z = 1, we have Here Keeping the order ω and ω 2 terms indicated, we have .
We are interested in the behavior of f a/A (η a /z, µ 2 ) for (1 − z) 1, so we neglect the z dependence of I(η a , µ 2 , z), which starts at order (1 − z) 2 . Then This gives us what Sterman and Zeng [46] call the "single power approximation." Sterman and Zeng argue that the single power approximation is numerically quite accurate in practice.
With the single power approximation, the ratios of parton distributions in eq. (7.21) becomes Then the Mellin moment variable appears in the exponent of our expressions, so that we can compare to standard results that are written in this form.

Rapidity dependence
We now consider the cross section to produce a muon pair with squared momentum Q 2 and rapidity Y : Here the lower limits on η a and η b are the momentum fractions at the Born level: (9.10) We are interested in the threshold region for Q 2 and Y , by which we mean that x a and x b are close to 1 and f MS a/A (η a , µ 2 ) and f MS b/B (η b , µ 2 ) are fast decreasing functions for η a > x a and η b > x b , respectively. In the threshold region, the flavor structure simplifies. Parton a must be a quark and b an antiquark, or vice versa. The flavor structure is carried by a function σ 0 that appears in the Born cross section: In the threshold region, we can write the parton level cross section in eq. (9.9) in terms of a dimensionless and flavor independent coefficient function C as Now we change integration variables to z and y defined in eq. (9.13): The limits − 1 2 log(1/z) < y < 1 2 log(1/z) come from the requirement that real emissions have positive components along the directions of p a and p b , so that η a > x a and η b > x b . Separately, the arguments x a e −y / √ z and x b e y / √ z of the parton distribution functions must be less than 1 or else the parton distribution functions will vanish. The requirements on the arguments of the parton distribution functions also implies that z > Q 2 /s.
We can now apply the single power approximation for the parton distributions, giving This implies that the parton distribution factor (and thus also the Born cross section) is maximum at that value of Y such that (N b − N a ) = 0. For our purpose of comparing to standard results, let us choose Y such that the parton factor is close to its maximum. Then N b ≈ N a . Now look at the integration in eq. (9.17). For the kinematic regime in which threshold summation is needed, (N a + N b )/2 is large. Then in the z-integration, z near 1 dominates. That means that the range of y is small. Since (N b − N a ) is small (compared to (N a + N b )/2), the factor exp(y(N a − N b )) can be approximated by 1, as argued in in refs. [21,22]. This gives 19) where N = (N a + N b )/2 and C(α s (Q 2 ), z) = These manipulations have given us a function of one variable to work with instead of a function of two variables. The function of one variable, C(α s (Q 2 ), z) is the function that appears in the rapidity-integrated cross section dσ/dQ 2 and is well studied. Eq. (9.19) gives us where C is the Mellin transform of C:

The standard result
Now we need C. We use the standard result [1,3] as given in eq. (3.1) of ref. [46] with factorization scale µ 2 f = Q 2 . We use the first term only in the cusp anomalous dimension, set the hard scattering function to 1, and omit the function D, thus dropping terms that contribute non-leading logarithms of N : (9.23)

The comparison
How does this compare to our results? In the parton shower approach, there are Sudakov factors K for the shower interval between the hard scattering that produces the muon pair and the first real parton splitting, then for the interval between the first real splitting and the second, and so forth. The integrands in the exponent of K decrease with decreasing splitting scale. Therefore we assume for the purposes of making a comparison that the first splitting occurs at a sufficiently small scale that we can just set that scale to zero and ignore the Sudakov factors for further splittings. This gives The factor K a comes from the exponential of V a − S a for initial state radiation from parton "a." The factor K b is the same expression applied to parton "b." The parton distribution functions in eq. (9.24) are those appropriate for a Λ 2 -ordered shower. They are related to the MS parton distribution functions by factors Z a and Z b defined in eq. (4.15), so that This matches the form of eq. (9.21). We need to check whether the factor C given by eq. (9.23) is the same as the leading approximation to Z a K a Z b K b .
We worked out the leading approximation to Z a in eq. (4.17): ) . (9.26) The factor K a comes from the exponential of V a − S a for initial state radiation from parton "a," integrated from an upper scale µ 2 a (t 0 ) = Q 2 , to a lower scale µ 2 a (∞) = 0. We use eq. (7.21) for K a . The term in eq. (7.21) proportional to ∆ ak is absent from the first factor K a for the Drell-Yan process and the term proportional to P reg aâ (z) can be omitted because it is not large. We also omit the i π term. Thus in eq. (7.21) we include only the main term, proportional to 1/(1 − z), we take the argument of the parton distributions to by fixed at Q 2 , and we approximate the lower endpoint of the z-integration as 1 − µ 2 /Q 2 instead of 1/(1 + µ 2 /Q 2 ). This gives In eq. (9.27) it is useful to interchange the order of integrations and change variables from θ(µ 2 ⊥ > m 2 ⊥ (a)) .

(9.28)
In the product Z a K a , the integrands combine in a nice way to give ) .

(9.29)
Although an analysis of other shower ordering schemes is beyond the scope of this paper, it is worth noting that if we had used k 2 ⊥ ordering or angular ordering for the shower, then Z a and K a would be different from what we have with Λ 2 ordering, but Z a K a would be the same. Now, apply the single power approximation (9.8). This gives We get the same factor for parton "b", so This matches the standard result (9.23) with two small changes. First, we have used a factor λ R in the argument of α s . This does not affect the leading logarithms. Second, we have 1 − z N instead of 1 − z N −1 . These forms are equivalent for large N . We conclude that the form of threshold logarithm summation that arises naturally in a parton shower is equivalent to the traditional forms that one gets in a direct-QCD analysis [8,9,46] or in soft-collinear-effective-theory as in ref. [37]. The shower form is less precise in that it does not allow an analysis beyond the leading approximation. On the other hand, it applies immediately to many processes with no further analysis.

Numerical comparisons
In this section, we exhibit two numerical tests of the threshold summation presented in this paper. In the first test, we look at the inclusive Drell-Yan cross section, p+p → e + +e − +X at 13 TeV. In the second test we look at the one jet inclusive cross section in proton-proton collisions at 13 TeV. We compare cross sections dσ, differential in whatever variables we choose, calculated with threshold summation to the corresponding cross sections without threshold summation.
The full cross section including threshold factors is dσ(full). This includes all of the terms in eq. (7.21) for (V a (τ )−S a (τ )) except the iπ term, with the color matrices calculated using the leading color approximation. 8 It also includes a factor Z a (η a , µ 2 f )Z b (η b , µ 2 f ), as defined in eq. (4.15). This factor relates the Λ 2 -ordered parton distributions to the MS parton distributions. Here µ 2 f is the factorization scale, characteristic of the hard scattering. For the jet cross section, α s at the hard interaction is evaluated at µ r = µ f . The parton shower then starts at scale µ 2 (t) = µ 2 f as given in eqs. (2.3) and (2.4). We will sometimes find it of interest to exhibit a cross section dσ(no ∆) in which we omit the term proportional to ∆ ak in eq. (7.21).
We can turn off all of the threshold effects to obtain a standard parton shower cross section dσ(std.).
In all of these cross sections the default calculation in Deductor begins with a factor where these functions obey the first order evolution equation (4.1). In the calculations presented in this section, we modify the initial parton factor to use NLO parton distributions by multiplying by a factor . The NLO parton distribution functions used are from the central CT14 NLO fit [77]. The parton distributions that obey the first order evolution equation (4.1) are simply obtained by using eq. (4.1) with the same starting distributions at the starting scale µ 2 start . For the Drell-Yan cross section, we find that R pdf is within about 5% of 1.
The Deductor (full) results depend on the parameter ψ min introduced in section 7.2.2. We use ψ min = 0.01. We have checked that varying ψ min by a factor 2 or 1/2 affects the cross sections examined by ±2% or somewhat less, depending on the observable.

Drell-Yan
In figure 1 we look at the inclusive Drell-Yan cross section, p + p → e + + e − + X at √ s = 13 TeV, as a function of the mass Q of the e + e − pair. We take µ f = Q. We show two curves from Deductor, one, dσ(full)/dQ, with the threshold effects turned on, the other, dσ(std.)/dQ, with the threshold effects turned off. Since the parton shower does not change Q, dσ(std.)/dQ equals the leading order (LO) perturbative cross section. For comparison, we show a perturbative next-to-leading order (NLO) calculation obtained from MCFM [78] with µ f = Q and CT14 NLO parton distributions. 9 We note that dσ/dQ decreases approximately exponentially as Q increases in the threshold region Q > 1 TeV. This reflects the fast decrease of the parton distribution functions as the momentum fraction increases. All three computed cross sections display the same approximately exponential behavior. However, the threshold correction has an effect that is large enough to notice even in this multi-decade semilog plot.
There is a theoretical uncertainty associated with the parton shower calculation, which we can estimate by changing the factorization scale µ f at which the initial parton distributions are evaluated and at which shower evolution starts. It is rather standard for the Drell-Yan cross section to choose the factorization scale to be µ f = Q. However, the maximum value of the transverse momentum of the e − or e + is Q/2, so, by analogy with jet production, for which µ f = P T (jet) is a widely used choice, µ f = Q/2 might seem a sensible choice here. On the other hand, one could choose µ f = 2Q. In figure 2, we plot the ratios 0.  of dσ(full) with these two scale choices to dσ(full) with µ f = Q. Based on this result, one might estimate a ±5% uncertainty. The precision of the Deductor calculation could be improved by matching to a NLO calculation of the Drell-Yan cross section, but we have not done this.
It is difficult to see small effects in semilog plots like figure 1, so, in figure 3, we show ratios K of cross sections to the Born cross section, dσ(std.)/dQ = dσ(LO)/dQ. In each case shown, we use µ f = Q.
We show first, in red, K(full), corresponding to the cross section with the threshold correction, dσ(full)/dQ. We see that the threshold correction is quite substantial and increases with Q.
We next show, in green, K(no ∆) for the cross section including just the part of the threshold correction in which we omit the term proportional to ∆ ak in eq. (7.21). The ratio K ∆ ≡ [dσ(full)/dQ]/[σ(no ∆)/dQ] is of some interest. For the Drell-Yan process, the ∆ ak term does not occur in the Sudakov factor between the hard interaction and the first real parton emission. After the first emission, there is color flow transverse to the beam so that the pattern of virtual gluon exchange is changed and ∆ ak can be non-zero. Since ∆ ak is itself proportional to α s , the perturbative expansion of K ∆ − 1 begins at order α 2 s . Thus we expect K ∆ to be close to 1. In fact, we find that K ∆ ≈ 1.03 for Q > 1 TeV. The parameter ψ min described in section 7.2.2 controls the integration range over which ∆ ak operates. As noted earlier, the cross section is sensitive to a factor 2 or 1/2 change in ψ min at a level of about ±2%, so the numerical value of K ∆ is not highly significant. What is significant is that K ∆ is indeed close to 1.
We display next, as a dashed black curve, the ratio K(NLO) of the NLO cross section to the Born cross section, as given by MCFM [78]. We note that K(NLO) is an increasing function of Q, as we might have expected since it includes some of the effect of threshold logs. We note also that the slope of the NLO curve remains rather constant as Q increases, in contrast to K(full), which has an increasing slope as the threshold logs build up.
Finally, we show as a purple, dashed curve, K obtained with the analytic threshold summation of Becher, Neubert, and Xu [37]. This curve is obtained by adapting the code for figure 8 of ref. [37] to CT14 parton distributions and √ s = 13 TeV. We regard the analytic B.N.X. curve as more precise than the Deductor (full) curve since the analytic result contains a high order of approximation in the summation of logarithms, while Deductor (full) is based on only a leading order parton shower. We note that, for Q > 2 TeV, the BNX curve agrees within about 3% with the Deductor result for dσ(full)/dQ.

Drell-Yan transverse momentum
In the previous subsection, we examined the Q dependence of the Drell-Yan cross section dσ/dQ, looking for the effects of steeply falling parton distribution functions when Q is large. Parton shower event generators can also predict the distribution of the transverse momentum Q T of the e + e − pair in the region of small Q T /Q, where logarithms of Q T /Q need to be summed. We have found [66] that a parton shower with virtuality based ordering, like Deductor, gives the same result at the next-to-leading-log level for logarithms of Q T /Q (without threshold logs) as the analytical treatment of ref. [79]. Now, with threshold effects included in a parton shower, we can examine both the logs of Q T /Q and the threshold effect at the same time, as in the analytical treatments of refs. [18,19,43,44]. We do not, however, have analytical knowledge of the level of accuracy of the parton shower treatment.
To study the Q T distribution at large Q, we examine We we divide by the integral of this over the Q T range 0 < Q T < 100 GeV to produce a distribution dN/dQ T normalized to  The normalized Drell-Yan transverse momentum distribution, dN/dQ T = (1/σ) dσ/dQ T for the LHC at 13 TeV. Here Q T is the transverse momentum of the e + e − pair, dσ/dQ T is dσ/(dQ T dQ) integrated over 2 TeV < Q < 2.1 TeV and σ is this cross section integrated over 0 < Q T < 100 GeV, so that the area under the curve is 1. The red curve that is lowest at small Q T is the full result with threshold effects, dN (full)/dQ T . The blue curve that is very slightly higher at small Q T is the result with no threshold effects, dN (std.)/dQ T . In these calculations, we have chosen µ f = Q. We also show the corresponding result obtained using ResBos [80,81] as a dashed black curve.
with the threshold correction omitted. We see that the threshold correction has only a very small effect. We also show the result dN (ResBos)/dQ T obtained with the analytical summation of logs of Q T /Q contained in ResBos [80,81]. The ResBos result does not contain a summation of threshold logs and so should be compared to Deductor (std.). The ResBos calculation contains smearing with non-perturbative functions that are fit to data. This smearing has not been included in Deductor. Thus it is not surprising that the Deductor distributions are somewhat narrower than the distribution from ResBos.
In figure 5, we examine directly dσ/dQ T defined in eq. (10.2) so that we can see the effect of the threshold factors on the normalization of the cross section. We take µ f = Q. We examine ratios K obtained by dividing dσ/dQ T by the cross section obtained with no threshold corrections, dσ(std.)/dQ T . We show two curves. In the upper curve, the numerator of K is the result with the full threshold correction, dσ(full)/dQ T . We see that there is a substantial, about 25%, threshold enhancement. This enhancement is weakly Q T dependent, 10 increasing from 23% to 26% over the range 0 < Q T < 100 GeV. We examine where this Q T dependence comes from by plotting also the ratio K obtained using dσ(no ∆)/dQ T in the numerator. Without the ∆ contribution in the Sudakov exponent, the threshold enhancement is flat as a function of Q T . Thus the small Q T dependence seen in dσ(full)/dQ T comes mainly from the ∆ term in the Sudakov exponent. This is easy to understand. The ∆ term appears only after we have an initial state emission. Having an initial state emission gives a transverse momentum recoil to the e + e − pair, so larger Q T should have a positive correlation with a larger threshold factor from the ∆ term. We can offer two observations. First, the Deductor (std.) curve for dN /(dQ T ) 10 There is a strong QT dependence at about QT = 1 GeV. This arises from the minimum pT allowed for emissions in the shower and is not really physical. dσ/dP T [nb/GeV] One jet inclusive cross section Figure 6. One jet inclusive cross section dσ/dP T for the production of a jet with transverse momentum P T and rapidity in the range −2 < y < 2. The cross section is for the LHC at 13 TeV. We use the anti-k T algorithm [82] with R = 0.4. The lower, blue curve is dσ(std.)/dP T , obtained with no threshold effects. The red curve is dσ(full)/dP T , obtained with threshold effects. The dashed, black curve is an NLO calculation [84]. In each case, we take the renormalization and factorization scales and the starting scale of the shower to be µ f = µ r = P T . agrees nicely with the ResBos curve, considering that there should be differences from non-perturbative smearing. Second, the effect of threshold logs reflected in the Deductor (full) curves in figures 4 and 5 is small and its sign appears to us to be quite sensible.

Jets
We now examine the one jet inclusive cross section dσ/dP T in proton-proton collisions at √ s = 13 TeV as a function of the jet transverse momentum, P T , integrated over the rapidity range −2 < y < 2. The jet is defined using the anti-k T algorithm [82] with R = 0.4 with the aid of FastJet [83]. Notice that the cross section dσ(std.)/dP T obtained with a standard shower with no threshold corrections is not the same as the Born cross section dσ(LO)/dP T because partons generated in the shower from initial state splittings can become part of the jet, while partons generated as daughters of the starting final state partons can escape from the jet.
In figure 6, we display three versions of dσ/dP T as functions of P T . In each case, we take the renormalization and factorization scales and the starting scale of the shower to be µ f = µ r = P T . The lower, blue curve is dσ(std.)/dP T , obtained with the parton shower with threshold effects omitted. The solid red curve is dσ(full)/dP T , obtained with threshold effects. We see that the threshold effect is large enough that it is evident even in this semilog plot. The black, dashed curve is the result of a purely perturbative next-toleading order (NLO) calculation [84]. We note that the parton shower calculation including the threshold effect is quite close to the NLO result.
There is a fairly substantial theoretical uncertainty associated with the parton shower calculation. To estimate this uncertainty, we examine the effect of changing the scale µ f = µ r at which the initial parton distributions and strong coupling are evaluated and at which shower evolution starts. In figure 6, we used µ f = µ r = P T . However, the minimum value of the dijet mass in the Born process is Q = 2P T . Thus µ f = µ r = 2P T might seem a sensible choice. One the other hand, jet cross sections are sometimes evaluated with µ f = µ r = P T /2, so P T /2 might seem a sensible choice. In figure 7, we plot the ratios of dσ(full)/dP T with µ f = µ r = 2P T and with µ f = µ r = P T /2 to dσ(full)/dP T with µ f = µ r = P T . Based on this result, we estimate a ±30% uncertainty in dσ(full)/dP T . This uncertainty could be reduced by performing a showered calculation matched to the NLO calculation.
In figure 8, we turn to several calculations of dσ/dP T presented as ratios K to the perturbative Born cross section, dσ(LO)/dP T . In this figure, all cross sections are evaluated at µ f = µ r = P T .
The lowest, blue curve is K(std.), obtained using dσ(std.)/dP T , in which there is a standard shower but the threshold effects are turned off. We see that dσ(std.)/dP T is only about 60% of the Born cross section. Since the cross section is so steeply falling as a function of P T , just a small amount of P T leakage out of the jet because of showering makes the cross section substantially smaller.
We now include the threshold correction, plotting the ratio K(full). This gives the red curve. We see that the threshold effect is very large and multiplies dσ(std.)/dP T by a factor between 1.3 and 2 for P T > 1 TeV. This produces a result dσ(full)/dP T that ranges from 90% to 120% of the Born cross section for P T > 1 TeV. We also show, as a solid green curve, the factor K(no ∆) obtained by omitting the term proportional to ∆ ak in the Sudakov exponent. We see that, for the jet cross section, this term makes a contribution to the exponent that is not negligible.
We show also as a black, dashed curve, K(NLO) corresponding to the perturbative NLO calculation from figure 6.
There are analytic summations of threshold logarithms [8,10,17,23,28]. We have used the computer programs of Kidonakis and Owens [17] and of de Florian, Hinderer, Mukherjee, Ringer and Vogelsang [28] to produce the cross sections dσ(K.O.)/dP T and dσ(FHMRV)/dP T , respectively. In these programs, the all-order threshold effect is expanded to order α 2 s to produce the calculated cross sections. In the Kidonakis-Owens formulation of threshold summation, there is no dependence on the algorithm used to define the jet. The FHMRV calculation is more sophisticated and includes the dependence  on the jet algorithm (for which we use the anti-k T algorithm with R = 0.4). In figure 8, we plot K(K.O.) as the green, dashed curve and K(FHMRV) as the purple, dashed curve.
We note that the Deductor (full) curve is roughly 30% below the analytic FHMRV curve for P T > 2 TeV. We also note that the scale variation test in figure 7 suggests that the Deductor (full) curve in figure 8 should be regarded as being uncertain to ±30%. In fact, we see in figure 7 that changing the scales to µ f = µ r = P T /2 makes the Deductor (full) cross section roughly 30% bigger. Thus the level of agreement between the Deductor (full) and FHMRV curves seems not unreasonable. It would, of course, be desirable to improve the precision of the shower cross section. This can be achieved by matching the shower calculation to the perturbative NLO correction to the cross section, but we have not yet undertaken this task.
We can draw a further conclusion from these comparisons. The parton shower has a hard job to perform, since it needs to include two large effects that act in opposite directions: loss of P T from the jet from showering and also the threshold enhancement. It Ratios of one jet inclusive cross section Figure 8. Illustrations of the one jet inclusive cross section as in figure 6. We show the results of calculations of dσ/dP T by plotting ratios K of dσ/dP T to the perturbative Born cross section dσ(LO)/dP T . The scales in all cross sections are µ f = µ r = P T . Going from the lowest to the highest curves at P T = 2 TeV, the lowest curve is K(std.) using in the numerator the showered cross section obtained with threshold effects turned off. The next is K(no ∆) obtained with the ∆ ak contribution turned off. The next curve is K(full) using the showered cross section with threshold effects. The next is K(LO) = 1. The next curve is K(NLO) using the perturbative NLO cross section. The next highest curve is K(K.O.) using the cross section obtained using the Kidonakis-Owens threshold effects program [17]. The highest curve is K(FHMRV) using the cross section calculated with the de Florian, Hinderer, Mukherjee, Ringer and Vogelsang algorithm [28]. seems to us remarkable that the calculation works to within the expected uncertainty.

Summary of the analysis
We have viewed parton shower evolution as the solution of equation (3.1), where ρ(t) represents the probability distribution of parton flavors and momenta and the density matrix in the quantum color and spin space in a statistical ensemble of event generation trials [67][68][69][70][71][72][73][74]. In this paper, we have ignored spin but still consider a full treatment of quantum color, even though in an actual implementation as computer code one has to make some approximations with respect to color. The shower time t is the negative logarithm of the hardness scale µ 2 considered. The shower starts at a hard interaction and evolves to softer scales. At shower time t, interactions that are softer than µ 2 are regarded as unresolvable, so that it is not meaningful to measure properties of the states described by ρ(t) at a finer scale. In the evolution equation, H I (t) represents real parton splittings, while S(t) leaves the number of partons, their momenta, and their flavors unchanged. Both operators are order α s ; we do not examine higher order contributions.
The operator S(t) gives us the Sudakov factor exp(− t 2 t 1 dt S(t)) that appears between two parton splittings. We have argued that S(t) should consist of two parts, as in eq. (6.6), Here F(t) represents the parton distribution factor in ρ(t) , so that the second term in S(t) gives the effect of changing the scale parameter in the parton distributions. The first term, S pert (t), accounts for order α s graphs that leave the number of partons unchanged. That is, S pert (t) represents one loop virtual graphs. We have (approximately) calculated S pert (t) in this paper.
In order to construct a parton shower based on eq. (11.1), one can use a trick. One can replace S(t) by V(t), where V(t) is constructed from the splitting operator H I (t) in such a way that the Born level cross section contained in ρ(0) is exactly conserved by the shower evolution. In fact, parton shower algorithms are typically based on V(t) instead of S(t). The difference [V(t) − S(t)] corrects −V(t). If we approximate the color using the leading color (LC) approximation or the LC+ approximation [70], then the color matrices are diagonal and exp( gives us a numerical weight factor that adjusts the cross section. We found that S pert (t) contains factors ±iπ times certain color matrices. These terms are very well known. (See, for example ref. [85].) As noted in ref. [70], these terms conserve the Born level cross section and could be included in V(t) instead of S(t). Although the iπ terms are of considerable physical interest, they are of secondary interest in this paper and have not been included in our numerical results.
The integrand in [V(t) − S(t)] contains ratios of parton distribution functions. The most important term has the form integrated over z in a small range near z = 1. This term is large if f a/A (η a /z, µ 2 ) falls steeply as 1 − z increases. Furthermore, f a/A (η a /z, µ 2 ) does fall steeply as 1 − z increases when η a is large, say bigger than 0.1. Thus a straightforward analysis of parton shower evolution leads us to the conclusion that there can be large corrections to the Born level cross section for a hard process. These contributions are naturally summed in a parton shower algorithm that incorporates S(t) because [V(t)−S(t)] appears as part of the Sudakov exponent.
Of course, these contributions are not summed if the Sudakov factor is exp[− dt V(t)], as is customary in parton showers (including ours, Deductor 1.0). In this paper, we have presented results from Deductor 2.0 [75], in which [V(t) − S(t)] is included within the LC+ approximation. (In our numerical results, we used the LC approximation.) The effects that arise from the term in eq. (11.3) are clearly connected with the effects of what are usually called threshold logarithms, which have been extensively studied. It thus is puzzling that the formulation in this paper contains ratios of parton distribution functions in an exponent, whereas standard threshold summation results never, to our knowledge, contain such factors. How, then, can these formulations be connected?
The answer can be understood in the treatment of the DGLAP parton evolution equation in the formulation of a parton shower. In a parton shower describing hadron-hadron collisions, at hardness scale µ 2 , one needs to include "unresolvable" initial state interactions as parton distribution factors, f a/A (η a , µ 2 ) and f b/B (η b , µ 2 ). When we come to a softer scale, we need to cancel the previous parton distribution functions and supply new ones. For this purpose, we can use That is, making use of the first order evolution equation,

(11.5)
This is what appears in shower evolution algorithms [48]. On the other hand, if we use the Mellin transform of the parton distributions, eq. (9.1), we havẽ Here γ(N ) it the matrix obtained by taking the Mellin transform of the evolution kernels. There are no parton distribution functions in the exponent. It is this formulation, or variations on it that do not work directly with the Mellin transformed evolution kernels, that typically appears in threshold summation calculations. For general purposes, eq. (11.6) is more powerful than eq. (11.5) because one does not need to have already solved the evolution equation to use it. However, eq. (11.6) is not more, or less, accurate than eq. (11.5).
In fact, we found in section 9 that a standard threshold summation in the case of the Drell-Yan process matches what the analysis of this paper gives when we use a leading saddle point approximation to connect the two results.
We now comment on the role of parton distributions in the formalism that we have presented. The ρ(t) in eq. (6.1) contains a factor representing the parton distribution functions at the resolution scale corresponding to t. In eq. (6.2), we defined an alternative statistical state ρ pert (t) in which this parton distribution factor has been removed. Then ρ pert (t) obeys the evolution equation (6.4) in which the parton distribution functions do not appear. At the end of the shower, we obtain the ordinary statistical state ρ(t f ) by multiplying by parton distribution functions appropriate to a low k 2 ⊥ scale, m 2 ⊥ (start), at which the shower turns off.
With this formulation, the whole hard scattering plus parton shower is a single perturbative process for which the resolution scale is m 2 ⊥ (start). There are no parton distribution functions except at the low scale. Now, the actual code works with ρ(t) and does use parton distribution functions at the hard scattering and at each shower stage. However, this use of parton distribution functions is only a trick [48]. Actually, all of the parton distribution functions approximately cancel except for those at the final low scale. 11 In order to make this cancellation work, Deductor 2.0 matches parton evolution to shower evolution by using parton distribution functions f a/A (η a , µ 2 Λ ) that evolve from the starting scale according a modified leading order evolution equation (4.8). It would, of course, be better to use an appropriate next-to-leading order evolution equation for the parton distribution functions, but we do not have a next-to-leading order shower algorithm. Thus we are stuck with a parton shower summation of logarithms based on leading order perturbation theory.

Choices in Deductor
A user of Deductor 2.0 has some choices.
The default choice is the calculation described above. With this choice, the parton shower sums threshold logarithms as described in this paper.
Another choice would be to eliminate the summation of threshold logarithms. This is easy to do, as described in section 10.
There is a third possibility. The user may want to retain the threshold logarithms but modify the way parton distribution functions appear in the calculated cross section. The default result for a cross section has the form where Q 2 is the scale of the hard scattering and some of the summation of threshold logarithms is contained in the factor .

(12.2)
See section 9.4. Here the parton distribution functions f MS a/A (η a , Q 2 ) are obtained from the parton distribution functions at the scale m 2 ⊥ (start) using the first order evolution equation (4.1).
Suppose that the user is only slightly interested in the details of the final state that a parton shower naturally specifies. Rather, the user is most interested in the hard scattering that initiates the parton shower. This user wants to have a calculation of the inclusive hard scattering cross section, including threshold corrections, that is as accurate as possible. Such a user might not want to have a cross section based on parton distributions f MS a/A (η a , Q 2 ) and f MS b/B (η a , Q 2 ) at the hard scale, since these parton distributions have been obtained by lowest order evolution from m 2 ⊥ (start). Instead, this user might prefer parton distributions f MS,NLO (η a , Q 2 ) that have been obtained with next-toleading order evolution. That is easily arranged by multiplying the default dσ in eq. (12.1) by a suitable weight factor: .

(12.3)
We have, in fact, done this in our numerical comparisons in section 10.

Outlook
We have presented a formulation of parton shower event generators in which the "threshold" enhancements of the cross section at large hardness scale Q 2 are included within the parton shower. This has a disadvantage compared to analytical summations of threshold logs: as presented here, the calculation has not been systematically extended beyond the leading logarithm approximation, whereas many of the analytical results are for a much improved order of approximation. However, the parton shower formulation has the advantage compared to analytical calculations that the same algorithm works for a wide variety of physical observables. As long as the desired Born level hard process is included in the parton shower code, the user simply has to specify the observable that is to be measured at the end of the shower. Compared to standard parton shower formulations that do not include threshold effects, the methods presented here have the advantage that they make the parton shower more accurate in a base level approximation in which matching to an NLO calculation has not been applied.
Every parton shower program is a little bit different. We have presented the threshold algorithms as needed for our program, Deductor. However, we believe that the methods presented here can be adapted with not much difficulty to other parton shower event generators.
were at the Munich Institute for Astro-and Particle Physics program "Challenges, Innovations and Developments in Precision Calculations for the LHC." It was completed while one of us (DS) was at the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara, program "LHC Run II and the Precision Frontier" which was supported by the U. S. National Science Foundation under Grant No. NSF PHY11-25915. We thank the MIAPP and the KITP for providing stimulating research environments. We thank Thomas Becher, Timothy Cohen, Hannes Jung, and George Sterman for helpful conversations. We thank Frank Petriello for advice about the Drell-Yan cross section and its perturbative calculation. We thank Jeff Owens for providing us with the Kidonakis-Owens code for threshold corrections to the one jet inclusive cross section. We thank Werner Vogelsang for providing us with the de Florian, Hinderer, Mukherjee, Ringer, Vogalsang code for threshold corrections to the jet cross section.

A Notation and kinematics
In this appendix, we collect some notations used throughout the paper and put them in one place. We only sketch the notation that we need. Details can be found in refs. [67][68][69][70][71][72][73][74].
In the parton shower, the partons are described by a state vector ρ(t) that represents the probability distribution of parton flavors and momenta and the density matrix in the quantum color space. (Recall that in this paper we effectively ignore quantum spin by summing over spins of the daughter partons after each splitting and averaging over the mother parton spin.) In this paper, all of the partons except top quarks are massless. We expand ρ(t) in basis states {p, f, c , c} m . This basis state represents two initial state partons with labels "a" and "b" and m final state partons with labels l. The partons have momenta p, flavors f , and colors c , c. The momentum fractions of the initial state partons are η a and η b . In the context of a parton splitting, we generally use p l to denote the momentum of parton l before the splitting andp l to denote the momentum of parton l after the splitting.
When a final state parton labelled l splits into partons l and m + 1 with momentâ p l andp m+1 , we characterize the splitting by its virtuality (p l +p m+1 ) 2 . Similarly, when the initial state parton with label "a" splits in backward evolution to a new initial state parton "a" and a new final state parton with label m + 1, we characterize the splitting by its spacelike virtuality (p a −p m+1 ) 2 .
The shower time that we use is related to the virtuality of a splitting: Here Q 0 is a fixed vector equal to the total momentum of all of the final state partons just after the hard scattering that initiates the shower. With this notation, the virtuality in an initial state splitting is It is convenient to use a dimensionless virtuality variable y = (p l +p m+1 ) 2 2p l · Q final state, Here Q = p a + p b is the total momentum of the final state partons just before the splitting. Then y is related to shower time by Since 2p a · Q = Q 2 , we have a convenient identity for µ 2 a (t): We sometimes use a squared transverse momentum variable k 2 ⊥ for a splitting of initial state parton "a." With the exact kinematics and momentum mappings used in Deductor, the emitted parton has label m + 1 and momentump m+1 . The part ofp m+1 orthogonal to the momenta of both incoming partons after the splitting,p a andp b , isp ⊥ m+1 , whose square is Here z is defined by z = η a /η a , where η a is the momentum fraction of parton "a" before the splitting (in backward evolution) andη a is its momentum fraction after the splitting. We note that the condition (p ⊥ m+1 ) 2 ≥ 0 implies The quantity (p ⊥ m+1 ) 2 vanishes when the emitted parton is collinear to parton "a," which corresponds to y → 0 with fixed (1 − z). However, it also vanishes when the emitted parton is collinear to parton "b," which corresponds to (1 − z − zy) → 0 with fixed y. For our purposes, we prefer a variable that matches (p ⊥ m+1 ) 2 in the collinear limit, but does not vanish in the anticollinear limit (1 − z − zy) → 0. We thus define a transverse momentum variable This discussion of k 2 ⊥ may be compared to that in section 2.3.1 of ref. [56]. We will need some standard flavor dependent constants. The number of flavors is N f . We have color factors C A = N c (with N c = 3) and C F = (N 2 c − 1)/(2N c ). We use f for a parton flavor, f ∈ {g, u,ū, d, . . . }. We define C f by In the kernels for evolution of parton distributions, constants γ f appear, with where T R = 1/2. We will use three momenta p l in a reference frame in which the total momentum Q of the final state partons before the splitting has zero space components. (Thus Q = p a + p b .) In the case of an on-shell virtual particle, p 0 l = | p l |. We then write with v 2 l = 1. We use E Q for the energy component of Q: Q = (E Q , 0 ). We often use the convenient shorthand In the case of an initial state parton, this is a a = a b = 1. We also sometimes use the shorthand ψ kl = 1 − cos θ kl where θ kl is the angle between a pair of partons with labels k and l.

B Calculation of the probability conserving integrand V(t)
For each real splitting in a parton shower, we need a Sudakov factor associated with the evolution of the system from the time t 1 of the previous splitting and the time t 2 of the new splitting. The standard form for this factor is T exp[− , where V(t) is the probability per unit dt to have a splitting at time t. Thus the Sudakov factor is the probability not to have had a splitting between t 1 and t 2 . This is the structure of the Sudakov factor used in Deductor 1.0.
In this paper, we calculate a numerical factor , where S(t) gives the effect of parton evolution and virtual graphs. We are not able to calculate this factor exactly. Rather, we calculate V(t) and S(t) for large t, or small y ∝ e −t as defined in eq. (A.4). In this appendix, we calculate V(t) for small y. Fortunately, although V(t) in Deductor is quite complicated, it has a simple structure for small y. Our aim in this appendix is to exhibit enough details of the calculation to enable the reader to reproduce it.
The operator V(t) is a sum over contributions from the partons in existence after time t 1 , as in eq. (5.4) [67]: (B.1) We will first examine the case of final state partons, with labels l = 1, . . . m. Then we will turn to the initial state partons, with labels "a" and "b."

B.1 Final state partons
We write V l (t) as [67] The sum includes all parton labels k = a, b, 1, . . . , m. When the operators V lk (t) act on a partonic state {p, f, c , c} m , they have the structure × {p, f, c , c} m .
As explained in section 5, the operators [(T l · T k ) ⊗ 1] act on the space of color density matrices for the parton state [67,70], with basis elements {c} m {c } m . A color generator matrix T a acting on parton l in the ket state multiplies a generator matrix T a acting on parton k in the ket state. The dot product indicates a sum over a = 1, . . . , N 2 c − 1. In [1 ⊗ (T l · T k )], we have the same construction with the color generators acting on the bra state.
From eq. (5.28) of ref. [70], we find that the functions λ V lk has the structure

(B.4)
Consider the first line on the right hand side of eq. (B.4). There is an integration over the variables that define the splitting of a mother parton with momentum p l into daughter partons with momentap l andp m+1 . In this paper, we take all partons to be massless. We will use the auxiliary variables We define the momentum fraction in the splitting bŷ where the auxiliary lightlike vectorñ l is We also define an azimuthal angle φ of the splitting using the part k ⊥ ofp l that is orthogonal to p l andñ l . The choice of flavors of the daughters can be specified by giving the flavor f m+1 of daughter parton m + 1. This gives us splitting variables y, z, φ,f m+1 . We can write the daughter parton momenta in terms of y, z, φ usinĝ The magnitude of the transverse momentum k ⊥ is given by Using eq. (8.20) of ref. [67], we find that integration over the splitting variables between shower times corresponding to y values y 1 and y 2 is accomplished with (B.11) The delta function that specifies the shower time is, from eqs. (A.1) and (A.3), (B.13) The first line inside the square brackets here is for a gluon self-energy graph with a quark loop. The color factor is T R = 1/2, but we follow the notation of eq. (B.3), in which this term comes with a color operator (T l · T l ) ⊗ 1 = C A 1 ⊗ 1 or 1 ⊗ (T l · T l ) = C A 1 ⊗ 1. The factor C A does not really belong here, so we remove it by dividing w ll by C A . The second line of eq. (B.4) covers gluon emission in a cut self-energy graph, while the third line covers gluon exchange between two lines, l and k.

(B.19)
Here x = λ(y) 1 + y z + 2a l y (1 + y)(1 + y + λ(y)) , . (B.20) Here we have changed conventions compared to ref. [68] and exchanged z ↔ (1 − z) and x ↔ (1 − x). This expression is rather complicated, but we only need its y → 0 limit: We also need w ll ({p,f } m+1 ) − w eikonal ll ({p,f } m+1 ) for g → g + g splittings,f l = g. It is given in eq. (2.50) of ref. [68]: We need only the y → 0 limit of this: Now we need the interference terms. Using eq. (5.3) of ref. [70] for w lk and eq. (7.12) of ref. [69] for the partitioning function A lk , we have To evaluate this, we write the spectator momentum after the splitting,p k , aŝ where ξ is a boost angle that is related to the angle θ lk between p l and p k in the Q = 0 frame by The parameter ω(y) is an additional boost angle, and in the case that k = a or k = b, ω(y) = 0. Since ω(y) is of order y and we are interested in the limit of small y, we will replace ω(y) → 0. The normalization parameter A k drops out of our calculation. The product in eq. (B.24) is quite complicated in general. However, it is reasonably simple in the limit of small y. We find Here we have noted that the terms proportional to y are negligible for small y unless (1−z) is small. Therefore in the coefficients of y and y(1 − z), we have replaced z by 1. We can integrate this over φ with the result [(1 − z) + a l y] 2 + 4a 2 l y 2 e 2ξ (1 + e 2ξ ) 1/2 .

(B.29)
We can then integrate this over z and take the small y limit of the result. We get We insert this result back into eq. (B.13). The integrals over φ and z of the other terms are simple, giving in the small y limit, We can perform the sum over flavorsf m+1 . The first term applies only if f l = g. When f l = g, there are N f equal terms. The third term applies if f l = g and there is one term. The second term applies when f l = g and there is one term. The fourth term applies for any f l . Thus (B.32) Using the constants γ f and C f from eqs. (A.10) and eq. (A.9), this is (B.33) When we insert this result into eq. (B.3), we have × {p, f, c , c} m .

(B.34)
In the second line of the right hand side, we can use k =l T l · T k = −T l · T l . Then in lines 1 and 2 we can use T l · T l = C f l . This gives × {p, f, c , c} m .

B.2 Initial state partons
We write V a (t) as The sum includes all parton labels k = a, b, 1, . . . , m. When the operators V ak (t) act on a partonic state {p, f, c , c} m , they have the structure as in eq. (B.3). We find from eq. (5.28) of ref. [70], In an initial state splitting, parton "a" with momentum p a becomes a new initial state parton with momentump a and a final state parton with momentump m+1 . The momentum of the other initial state parton is unchanged: where n ⊥ ·p a = n ⊥ ·p b = 0 and n 2 ⊥ = −1. Note that y =p a ·p m+1 /p a ·Q, where Q = p a +p b , as in eq. (A.3). This kinematics requires y < (1 − z)/z, or Eq. (B.38) contains a ratio of parton distributions, which are evaluated at momentum fractions defined by p a = η a p A andp a =η a p A . We take the scales of the parton distributions to be the virtuality µ 2 a (t) defined in eq. (A.2). There are also factors n c (a) that count the number of colors (3 or 8) carried by partons of flavor a.
Using eq. (8.20) of ref. [67] together with eqs. (A.28) of ref. [73], we find that integration over the splitting variables between shower times corresponding to y values y 1 and y 2 is accomplished with (B.41) The delta function that specifies the shower time is, from eqs. (A.1) and (A.3), There are four terms here. The first three are for direct graphs. The first is for (â, a,f m+1 ) = (g, q,q) and (â, a,f m+1 ) = (g,q, q). The second is for (â, a,f m+1 ) = (q, g, q) and (â, a,f m+1 ) = (q, g,q). Here the color factor is T R = 1/2, but this term multiplies T a · T a = C A , so we need to divide by C A . The third term is for (â, a,f m+1 ) = (g, g, g), (â, a,f m+1 ) = (q,q, g) and (â, a,f m+1 ) = (q, q, g). Here there is a soft gluon singularity, which is subtracted. The fourth term is for interference graphs, in which a gluon is exchanged between parton "a" and parton k.
Now we need the splitting functions. We take the limit y 1. From ref. [73], eq. (A.42), we have for (â, a,f m+1 ) = (g, q,q) and (â, a,f m+1 ) = (g,q, q), From ref. [73], eq. (A.45), we have for (â, a,f m+1 ) = (q, g, q) and (â, a,f m+1 ) = (q, g,q), From ref. [73], eqs. (A.35) and (A.39), we have for (â, a,f m+1 ) = (g, g, g), (â, a,f m+1 ) = (q,q, g) and (â, a,f m+1 ) = (q, q, g), The eikonal function is defined in eq. (2.10) of ref. [68] as Using eq. (B.39) gives The term z/(1 − z) here removes the 1/(1 − z) singularity from P aa (z). The third term, proportional to y has a 1/(1 − z) singularity. This can give a log(y) contribution to an integration over z down to (1 − z) = y. However, we can neglect y log(y). Thus we can throw this term away. This gives Now we need the interference terms. From ref. [70], eq. (5.3), we have This multiplies A ak . We use eq. (7.12) of ref. [69]: The product is To proceed further, we note that we need the function in eq. (B.54) in the limit of small y. The momentump k is related to the momentum p k by a Lorentz transformation that becomes the unit operator when y → 0. Thus we can neglect the difference between p k , which varies as we integrate over splitting variables (y, z, φ), and p k , which is fixed. For this reason, we substitute p k forp k in eq. (B.54). Then we use Here A k is a normalization factor that cancels in eq. (B.54) and θ ak is the angle between the three-vector parts of p k and p a in the frame in which Q has only a time component, while u ⊥ is a vector transverse to p a and p b with u 2 ⊥ = −1. Using the parameterizations in eqs. (B.39) and (B.55), we find (1 − cos θ ak ) 2(1 + cos θ ak )yz 2 + (1 − cos θ ak )(1 − z) − 2z sin θ ak cos φ y(1 − z − yz) .

(B.56)
We can perform the averaging of this over φ exactly: (B.57) We can write this in a suggestive form as We are interested in this in the small y limit, in which y 2 z 2 /Ψ ak (z) 2 1. Then the second term in the denominator of eq. (B.58) is non-negligible only when z is close to 1. Thus we can replace Ψ ak (z) by ψ ak defined in eq. (A.13): We may note that the angle θ ak is small when partons "a" and k are the daughter partons of a previous initial state splitting that was nearly collinear. It is allowed in Deductor to have an initial state splitting with small momentum fraction z ak , so that the new initial state parton "a", which is the mother parton for the next splitting, has a much larger momentum fraction than the previous initial state parton. In this case, the virtualities of the initial state partons in successive splittings are not strongly ordered, so that one can have y > 1 − cos θ ak . This regime of multi-regge kinematics is discussed in section 5.4 of ref. [72]. This is the opposite kinematic regime from that of threshold logarithms, so we ignore this possibility in this paper. However, we still use y 2 z 2 instead of just y 2 in the denominator of eq. (B.61) in order to keep the result reasonably accurate even when y 1 − cos θ ak .

(B.63)
In the first term in eq. (B.63), we can replace (T a · T a ) by C a . We divide the second term into three terms by defining ∆ ak (z, y) according to 1 × {p, f, c , c} m . We note that for k = b, we have cos θ ab = −1 and 1/ψ 2 ab = 0. Then the integral is just − log(y). This gives × {p, f, c , c} m .

(B.67)
Now there are three k = a terms. The first has the good feature that its only k dependence is in the color factor, so that we can sum it over k. In the second term, we note that by its construction, ∆ ak (z, y) has a 1/(1 − z) singularity for (1 − z) y/ψ ak but is suppressed compared to 1/(1 − z) when (1 − z) y/ψ ak . Thus the second term has the good feature that, because of the structure of ∆ ak (z, y), the only important integration region for the z-integration is 0 < (1 − z) y/ψ ak . We also note that 1/ψ ak = 0 when cos θ ak = −1. In this limit, ∆ ak (z, y) = 0. We have cos θ ak = −1 for k = b, so ∆ ab (z, y) = 0. This eliminates one term in our sum over k. The third term has the good feature that we have been able to integrate it, at least approximately.
In the second term in eq. (B.67), we can sum over k using k =a (T a ·T k ) = −(T a ·T a ) → −C a . In the last term, we can separate the log(1/y) contribution and perform the color sum in the same way. This gives × {p, f, c , c} m . (B.68) We will use this result in eq. (5.6).

C Calculation of the virtual graphs
The Sudakov exponent S that appears in eq. (6.1) consists of two terms, as given in eq. (6.6): a term S pert (t) that comes from virtual graphs and a term that accounts for the evolution of the parton distribution functions. In this appendix, we outline the calculation of S pert (t).
Consider first an operator S pert tot that corresponds to the one loop virtual graphs that contribute to shower evolution. (More precisely, because of the minus sign in eq. (6.1), S pert tot corresponds to the negative of the one loop virtual graphs.) In these graphs, we integrate over a loop momentum k. There will be ultraviolet divergences that come from one loop corrections to QCD propagators and vertices. We suppose that these are removed by renormalization. There may be an additional ultraviolet divergence that arises from letting the scale of the loop momentum be much larger than the scale Q 2 of the hard interaction that initiates the shower. This can happen if we make the approximation Q 2 → ∞ inside the graph. We should simply arrange to have a regulator in the integrand that eliminates such a divergence.
The operator S pert tot will also have infrared divergences, corresponding to the integration regions in which k → 0 or k becomes collinear with the momentum of one of the external lines of the graph. We regulate the infrared divergences with dimensional regulation or some other method.
We will represent S pert tot as an integral over a shower time t, as we have done for real emission diagrams: The shower time t corresponds to the negative of the logarithm of the hardness scale of the integrand, in analogy with the definition of t in real emission graphs. Thus the infrared divergences are associated with t → ∞. This general idea does not, however, tell us exactly how to define t. One way to proceed would be to use dimensional regulation and subtract the infrared poles. Then the result would depend on a parameter µ 2 . Then we could identify e −t with µ 2 /Q 2 and S pert with the derivative of the graphs with respect to log(Q 2 /µ 2 ). However, in eq. (C.24) below, we will make a more direct identification based on the physical meaning of our version of the shower time.
Using the notation analogous to that of eq. (5.4) for real emission diagrams, we can write The part associated with splitting of parton l is S pert l (t). As in eq. (B.2) for real emission graphs, we write S pert l (t) as S pert 3) The first term describes self-energy interactions. In the second term, a virtual gluon is exchanged between parton l and parton k. The sum includes all parton labels k = a, b, 1, . . . , m except for k = l. The operators S lk (t) have the color structure As in appendix B, T a k inserts a color generator matrix T a on line k. The functions S R lk (t) are just the complex conjugates of the functions S L lk (t), so we need only to define and analyze the functions S L lk (t). We will begin with the case of interference diagrams: k = l. We start with the case that both l and k represent final state partons. Then we will look at the case in which both l and k represent initial state partons. Finally we let one of the partons be in the initial state while the other is in the final state. Once we have covered interference diagrams, we look at self-energy diagrams: k = l.

C.1 Final state interference diagrams
Let's work with S L lk ({p, f } m ; t) in the case that l and k are (different) final state partons in the ket amplitude. We look at the contribution to S L lk ({p, f } m ; t) from gluon exchange between two partons. As we do throughout this paper, we take all partons to be massless. The momenta of partons l and k are p l and p k . The gluon carries momentum q from line k to line l, so that, inside the loop, line l carries momentum p l − q and line k carries momentum p k + q.
In analyzing the virtual graphs, we follow as much as possible the treatment of the real graphs that we have used in refs. [67][68][69]. In particular, this means that we calculate the virtual graphs in Coulomb gauge.

C.1.1 Integral in Coulomb gauge
We start with the full interference graph. We will want to identify the shower time, but we have not done that, so we start with the integral over t of the functions that we ultimately want. Also, when a gluon is exchanged between partons l and k, it is undefined which is the primary emitter and which is playing only a helping role. That is, we want S L lk to represent an emission from parton l with parton k as helper, while S L kl will represent an emission from parton k with parton l as helper. We have not yet defined how the total graph is partitioned into these two parts, so we begin with S L lk + S L kl . Thus we start with the definition (including the minus sign in eq. (6.1), so that we write the negative of the usual Feynman diagram), (C.5) We work in Coulomb gauge in a reference frame in which Q = 0. The integration variables are defined by q = (E, q ) in this frame. The integral is both ultraviolet and infrared divergent. For technical reasons have inserted an ultraviolet cutoff | q | < M , where we will take M to be large. The infrared divergence is associated with large positive values of the shower time t. We could imagine that the infrared divergences are regulated, but once we select a fixed value of t, the regulation is not needed. For this reason, we do not specify a regulation method. In Coulomb gauge at a fixed shower time t, only the soft integration region for q is important. For this reason, it is appropriate to use the eikonal approximation (as in Deductor). We have applied the eikonal approximation in eq. (C.5). The numerator of the gluon propagator in Coulomb gauge is whereq = (0, q ). Now, the integrand in eq. (C.5) is complicated because of the numerator D(q) of the gluon propagator in Coulomb gauge. However, we can simplify it by writing with an analogous definition of S L kk ({p, f } m ; t; eikonal). This is the eikonal approximation to the self-energy graph for parton l in Coulomb gauge. Recall that we associate S L lk with emissions from parton l and S L kl with emissions from parton k.
Since q · P lk = 0, none of the q dependent terms in D(q) contribute and we are left with − P lk · D(q) · P lk = P 2 lk = −2q · p l q · p k p l · p k . (C.12) Thus That is, we get the familiar dipole formula for one gluon exchange in Feynman gauge.

C.1.2 Dipole part
We now analyze the Feynman gauge eikonal integral in eq. (C.13). We write In the first term, the gluon propagates forward in time from parton k to parton l. In the second term, the gluon propagates forward in time from parton l to parton k. In the second term, we redefine q → −q, so that the direction of q is the direction of propagation forward in time. Then . .

(C.17)
We can immediately perform the E-integration to give .

(C.18)
Now, we can rewrite this using This gives That is Now consider the first term. It has essentially the structure of the graphs for the emission of a real gluon with q 2 = 0, although it is not quite the same as our real emission factors because it does not contain momentum mappings that allow an on-shell parton to split into two on-shell partons. The integrand has poles at | q| = q· v l and at | q| = q· v k . This reflects splittings both of parton l and of parton k. In order to separate these splittings, we multiply the integrand by 1 = A lk + A kl , where The term containing A lk is associated with S L lk ({p, f } m ; t), while the term containing A kl is associated with S L kl ({p, f } m ; t). In the iπ term, the two contributions are actually equal, but we can associate the first with S L lk ({p, f } m ; t) and the second with S L kl ({p, f } m ; t). This gives (C.23) We can now identify the shower time with This matches the definition (A.3) that we used for a real gluon emission. 12 We introduce this as a delta function, recognizing that dt · · · is equivalent to d log(y) · · · . Thus (C.25) We will want to apply different methods for the two integrations. 12 In the case of real gluon emission, we take q and l to be the parton momenta after the splitting. The momentum of parton l before the splitting is l + q, but here we use just l because q is small.

C.1.3 Dipole real part
We examine first the real part of the integral (C.25). We introduce transverse and longitudinal coordinates for q: where q ⊥ · p l = 0 and where we have defined as in eq. (A.12). We denote by φ the azimuthal angle of q ⊥ relative to the ( p k , p l ) plane.
We first need to find the value of | q ⊥ |. From eq. (C.24), we have This implies that the cutoff | q | < M amounts to Note that there is a minimum possible value of (1 − z), corresponding to q 2 ⊥ = 0: (1 − z) > a l y . (C.33) We can write the integration over q as To perform the q 2 ⊥ integration against the delta function, we note that This gives We can integrate this over φ: The log(M/| p l |) term will cancel against an identical term in the integral that we subtracted and have to add back.

C.1.4 Dipole imaginary part
Now we examine the imaginary part of the integral (C.25), which we rewrite slightly as ) .

(C.42)
We can immediately take the limit M → ∞. It will be useful to choose coordinates (ξ, η, λ) based on the orthogonal vectors ( v k + v l ) and ( v k − v l ): where u is a unit vector orthogonal to v k and v l . We then have d q = 1 1 − cos θ kl dξ dη dλ . (C.44) We can immediately perform the λ integration using δ( q · ( v k − v l )) = δ(λ) .

(C.55)
We can perform the E-integration to get We recognize the denominator as defining the shower time, so We put our contributions back together, inserting G L lk ({p, f } m ; t; dipole) from eq. (C.53) and G L ll ({p, f } m ; t; eikonal) from eq. (C.60) into eq. (C.7): Here the cutoff dependent terms proportional to log(M/| p l |) have cancelled. We will use this result in eq. (6.14).

C.2 Initial state interference diagram
Let's now look at the case that the active parton is one of the initial state partons, l = a, and the helper parton k is the other, k = b. Thus, we examine S L ab ({p, f } m ; t), corresponding to gluon exchange between the two initial state partons. The gluon carries momentum q from line "a" to line "b," so that, inside the loop, line "a" carries momentum p a − q and line "b" carries momentum p b + q.
We start with the eikonal approximation to the exchange in Coulomb gauge, 2p a · D(q) · p b (q · p a − i )(q · p b + i )(q 2 + i ) θ(| q | < M ) .

(C.62)
There is an ultraviolet cutoff | q | < M that we eventually remove. This integral is exactly the same as we had in eq. (C.5) for the final state case, with p l → p a and p k → p b . We can apply the same treatment, partitioning the integral into two terms using the partitioning function (C.22) and identifying the shower time using eq. (C.24). Thus we can simply use the result in eq. (C.61), noting that cos θ kl = −1: (C.63) We will use this result in eq. (6.15).

C.3 Initial state -final state interference
We now examine the case of gluon exchange between one of the initial state partons, say l = a, and a final state parton k. Thus, we examine S L ak ({p, f } m ; t) and S L ka ({p, f } m ; t). An exchanged gluon carries momentum q from line "a" to line k, so that, inside the loop, line "a" carries momentum p a − q and line k carries momentum p k − q.

(C.67)
As with eq. (C.15) we write this using E and q: (C.68) We can immediately perform the energy integration, noting that the term 1/(E+| q|−i ) does not contribute: (C. 69) Notice how this compares to the equivalent result for virtual gluon exchange between two final state partons, as in eq. (C.18), or two initial state partons. Here there is no imaginary part. The integrand has poles at | q| = q · v a and at | q| = q · v k . This reflects splittings both of parton "a" and of parton k. In order to separate these splittings, we multiply the integrand by 1 = A ak + A ka , where . (C.70) The term containing A ak is associated with G L ak ({p, f } m ; t), while the term containing A ka is associated with G L ka ({p, f } m ; t). Thus we define θ(| q | < M ) .

(C.71)
We can now identify the shower time with as we did for a final state splitting. We introduce this as a delta function, recognizing that dt · · · is equivalent to d log(y) · · · . Thus × 1 E Q y (| q| − q · v k + E Q y) .
(C. 73) where k ± are the (on-shell) parton momenta after the splitting. The virtuality (k + + k − ) 2 is (| k + | + | k − |) 2 − ( k + + k − ) 2 . This is normalized by dividing by 2p l · Q = 2p 0 l E Q . In a real emission, a small amount of momentum is taken from elsewhere in the event to put p l on shell. For the virtual graph, we simply replace p 0 l → | p l | = | k + + k − |. Thus we identify In the result from ref. [86], the integral is written as an integration over y, the azimuthal angle φ of k + around the direction of p l , and and a variable 13 The variableq 2 = 2| p l |E Q y appears instead of y in ref. [86]. We include here the factor 1/2 in eqs. (C.76) and (C.77). Ref. [86] gives directly the right hand side of these equations without the factor 1/2. Also, we maintain the notation of eq. (C.4) in which S L ll multiplies a factor [(T l · T l ) ⊗ 1]. For the gluon case, this is a factor C A [1 ⊗ 1] and for the quark case, it is a factor C F [1 ⊗ 1]. For this reason, we remove a factor C A or C F from the results as given in ref. [86]. As in previous sections, we will use the convenient abbreviation a l = E Q /(2| p l |) from eq. (A.12).
The variable y is proportional to e −t , where t is the shower time. Thus integrating over t is the same as integrating over y, with dt = d log y. We are interested in S L ll ({p, f } m ; t), the integrand of the log y integration. Integrations over x and φ will remain.
We now turn to the gluon and quark cases separately.

C.4.1 Gluon self-energy
For the gluon self energy, we find from ref. [86] S L ll ({p, f } m ; t; gluon) = The coefficients A T,J are (C.82) 13 In the quark case, k− is the momentum of the quark line inside the loop, so x is the momentum fraction of the gluon. However, we present the results symmetrized over x ↔ (1 − x), so that it does not matter whether x or (1 − x) is identified with the gluon in the loop. Ref. [86] explains how to remove the symmetrization, but we do not need to do that here.
Here we have renormalized the graphs according to the MS prescription, but have subtracted a numerical function that gives the same result as subtracting a pole. The parameter µ 2 R is the MS renormalization scale. We can perform the integrations. The φ-integral is trivial. Performing the x-integral gives (1 + 2a l y)(1 + 8a l y + 20a 2 l y 2 ) (1 + 4a l y) 3/2 log (1 + √ 1 + 4a l y) 2 4a l y .
(C. 84) We are interested in the small y limit of this. There is a constant term and a term proportional to log(y): (C.85)