Parton showers with more exact color evolution

Parton shower event generators typically approximate evolution of QCD color so that only contributions that are leading in the limit of an infinite number of colors are retained. Our parton shower generator, Deductor, has used an"LC+"approximation that is better, but still quite limited. In this paper, we introduce a new scheme for color in which the approximations can be systematically improved. That is, one can choose the theoretical accuracy level, but the accuracy level that is practical is limited by the computer resources available.


I. INTRODUCTION
Parton shower event generators for hadron collisions, such as Herwig [1], Pythia [2], and Sherpa [3], are essential for the analysis of experiments at the Large Hadron Collider. They treat QCD color in the leading color (LC) approximation, that is at leading order in an expansion in powers of 1/N 2 c , where N c = 3 is the number of colors. Previous versions of our parton shower generator, Deductor [4][5][6][7][8][9][10][11][12][13], use what we call the LC+ approximation [7], which includes some terms suppressed by powers of 1/N 2 c . This paper concerns improvements in the color treatment in Deductor that could be applicable to other parton shower generators. The treatment of parton spin is also important, but we simply ignore spin here. 1 It is also of interest to consider how color works in a parton shower at an arbitrary order of perturbation theory for the splitting functions that generate the shower [14]. However, we work only to first order in α s for the splitting functions here.
One can define an evolution equation for a parton shower with leading order splitting functions such that the evolution is exact in color. This is relatively straightforward [4], as we explain briefly below. Once we understand the exact evolution equation, we are faced with trying to implement it as a computer program. This is, so far as we know, impossible with any available computer. However, one should not be discouraged. What we really need is not a numerical answer that is exact with respect to color for the calculation of a cross section σ of interest, but rather an approximate answer that can be systematically improved.
What we seek is an algorithm for a parton shower that depends on parameters that control the level of approximation with respect to color. We can begin with the cross section calculated at the lowest level of approximation, call it σ(0). Then we can change parameters in the calculation so that we get successively better approximations, σ(1), σ(2), . . . . We will find that the successive calculations use more and more computer resources, so that there will be a practical limit to how exact we can be. Note that |σ(n)−σ(n−1)| serves as an error estimate for σ(n). If |σ(n) − σ(n − 1)| is not as small as we would like by the time that we run out of computer resources, we will have to admit the limitations of the calculation. However, even in that case, we obtain a calculation with an error estimate. With just σ(0), we have a less precise calculation with no error estimate.
There has been work on extending the accuracy of parton shower algorithms beyond the leading order in an expansion in powers of 1/N 2 c [7,11,[15][16][17]. To date, however, what we see as the main obstacle to a systematically improvable treatment of color has not been overcome. This obstacle has been nicely stated in Ref. [17]: "To fully include all subleading N c terms in the soft and collinear limits, virtual color rearranging terms associated with the same singularity structure should also be kept. To accomplish this, a full resummation of virtual exchanges is needed. Unfortunately, within the current event generator structure these contributions cannot be included, and we postpone their inclusion for future work." The aim of the present paper is to provide an event generator structure that does include these contributions in a systematically improvable fashion.
Before we turn to approximate versions of color in a parton shower, we briefly describe a shower evolution equation that is exact in color, as in Ref. [4]. We view shower evolution as an application of the renormalization group. At any stage, we include hard interactions and remain inclusive over softer interactions. Then the picture changes when we change the definition of where "hard" interactions end. We define t = log(µ 2 0 /µ 2 ), where µ 2 is a arXiv:1902.02105v1 [hep-ph] 6 Feb 2019 measure of hardness 2 in a splitting and µ 2 0 is a reference value of µ 2 . The shower develops with increasing t as the interactions that are included become softer and softer.
At any stage of the shower we have m final state partons plus two initial state partons, "a" and "b," with momenta and flavors {p, f } m = {p a , f a , p b , f b , p 1 , f 1 , . . . , p m , f m }. These partons carry QCD color, so we can describe them using a quantum amplitude M ({p, f } m ) , which is a vector in the color space of m final state partons plus two initial state partons. In order to keep the notation simple, we ignore spin in this paper, but if we were to include spin, then M ({p, f } m ) would be a vector in the color×spin space.
Shower evolution is probabilistic since, at any stage of the shower, we are integrating over probabilities for potential splittings that are too soft to be visible at that stage. For this reason, we use the language of quantum statistical mechanics to describe shower evolution. With this language, we do not use pure color amplitudes but rather a color density operator, We use the color density operator as the basis for parton shower evolution. We note that it is sometimes also used explicitly for analytical summations of large logarithms [18][19][20][21].
We can expand ρ({p, f } m ) in color basis states, Deductor uses the trace basis (which might better be called the color string basis) for color basis states {c} m [4]. The set of all functions ρ({p, f, c, c } m ) constitutes a vector space, which we call the statistical space. 3 We represent the function ρ as a ket vector, ρ . The rounded end of the ket is meant to distinguish a vector in the statistical space from a vector in the quantum color space. Note that ρ represents the whole function ρ, which gives the distribution of simulated events as a function of the number of partons and their momenta, flavors, and colors. A parton shower event generator generates particular events by Monte Carlo sampling from this distribution. In each event, there is, for instance, a certain choice for the number of partons and their momenta. 2 As reviewed in Ref. [13], the default hardness measure in Deductor is a variable Λ 2 that is proportional to the virtuality of the splitting. See Eq. (35). However, an alternative in Deductor is an appropriately defined transverse momentum squared k 2 T in the splitting. 3 To be more precise, ρ is a function of m and {p, f, c, c }m. Additionally, ρ depends parameters x derived from prior states, {p, f, c, c } m with m < m. For example, ρ depends on the vector Q 0 equal to the total momentum of the final state partons at the start of the shower, which is used to define the shower evolution variable. Normally we suppress the dependence on other parameters x in statistical states and in functions.
We find it useful to use a notation in which linear operators act on vectors ρ in the statistical space. Thus we might write Martínez, De Angelis, Forshaw, Plätzer, and Seymour [16] have recently analyzed the influence of color on parton evolution via soft gluon emission using the color density operator, but with a different notation from ours. If the operators A i in Eq. (3) have a suitable product form, these authors would write ρ and ρ as density operators as in Eq. (2) and write There is no physics difference between these two notations. We will use the notation of Eq. (3).

II. EVOLUTION EQUATION EXACT IN COLOR
We can now discuss the evolution equation that forms the basis for the approximations used in Deductor [4,5,7]. The statistical state changes as the hardness resolution varies, so that it is a function ρ(t) of t. We can write for any t and t 0 . There is also a threshold factor U V that appears at the start of the shower [13,14]. See Eq. (84). We return to this factor later. The shower evolution equation operator U(t, t 0 ) with full color obeys the evolution equation, The operator H I (t) creates a parton splitting, increasing the number of partons by one. The new parton carries color and the colors of the old partons change. The operator V(t) leaves the number of partons unchanged. It does, however, change the color state of the partons. It carries the color structure of virtual graphs. These two operators are related by an identity, Here, for any parton state ρ , to calculate 1 ρ , we integrate over all of the parton momenta in ρ , sum over flavors, and take the trace over colors. This gives us the total probability associated with ρ . Because of Eq. (7), together with the initial condition U(t 0 , t 0 ) = 1, we have for any statistical state ρ . This says that probability is conserved in the evolution of the state.
The operator V(t) contains a term that we can call V iπ (t) that we calculate from the imaginary part of virtual graphs. This operator has the form (assuming massless partons) given in Eq. (10.14) of Ref. [7], Here T a represents the insertion of a color matrix T c on incoming parton line "a", T b represents the insertion of a color matrix T c on incoming parton line "b", and the dot in (T a · T b ) represents a summation over the octet color index c. In [(T a · T b ) ⊗ 1] the color matrices act on the ket state, while in [1 ⊗ (T a · T b )], they act on the bra state. When we take the color trace, we get 1 V iπ (t) = 0. The operator V(t) also contains terms with real coefficients that reflect the color structure of the real parts of virtual graphs, Here [(T l ·T k )⊗1] inserts color matrices on the color lines of partons with indices l and k in the ket state, while [1 ⊗ (T l · T k )] inserts color matrices on the color lines of partons with indices l and k in the bra state. Here l or k or both can be the indices "a" and "b" of the initial state partons and k can equal l.

III. EVOLUTION IN THE LC+ APPROXIMATION
We can use the LC+ approximation described in Ref. [7]: Ref. [7] does not give an LC+ approximation for V iπ (t). 4 In this paper, we simply take 5 The LC+ operators H lc+ (t) and V lc+ (t) obey This gives us so that probability is conserved in LC+ evolution. 4 The suggestion in Ref. [7] that the first term in Eq. (10.8) of that paper might be included in the LC+ approximation does not work because Eq. (13) below would fail. 5 One can, at least in principle, include V iπ (t) in the LC+ approximation. We leave the exploration of this possibility to future work.
The differential equation Eq. (11) can be solved in the form Here N lc+ (t 2 , t 1 ) is the no-splitting operator, It is well to recall here an essential point: the operator V lc+ (τ ) is diagonal in the color basis that we use, the trace basis, so that it is practical to calculate its exponential.

IV. EVOLUTION BEYOND THE LC+ APPROXIMATION
Now, what if we want shower evolution with full color? Then we need where Note that since we have set V lc+ iπ (t) = 0, V iπ (t) is included in ∆V(t): This differential equation can be solved in the form Here N (t 2 , t 1 ) is the no-splitting operator, where T denotes time ordering for the non-commuting operators in the exponent, with latest times to the left. Unfortunately, we cannot include ∆V(τ ) here in exponentiated form because of its complicated color structure. Thus we write the evolution equation for N (t, t 0 ): We can solve this in the form This can be solved iteratively. The first three terms are It is convenient to write the solution as where When we use X (t, t 0 ), we understand that it is expanded to whatever order in ∆V(τ ) that we need.
Using the operator X , we write the evolution equation Eq. (20) in the form This generates a shower, at least in principle. At each step in the shower the splitting operator is Notice that X (τ, t 0 ) is an operator on the color space, but does not create any new partons. Once the color is transformed by X (τ, t 0 ), the operator H lc+ (τ )+∆H(τ ) creates a new parton and further modifies the color. The net operator in Eq. (28) is then a splitting operator in the sense that it creates a new parton. It is also a non-trivial operator on the color state.
Our goal now is to treat the operator ∆H(τ ) and the operator ∆V(τ ) in X perturbatively. To avoid confusion, we note that expanding in powers of ∆H(τ ) and ∆V(τ ) is not equivalent to expanding in powers of 1/N 2 c . In order to make the evolution equation (27) practical for a computer program, we will need to rearrange it. To do that, we first review the singularities that control the evolution and analyze the different roles of soft and collinear singularities.

V. SPLITTING VARIABLES AND SINGULARITIES
The splitting operator H I (t) is singular in the limit of very large shower times t, corresponding to very small splitting virtualities. In order to study this limit, it is convenient to define a dimensionless virtuality variable y. For a final state splitting in which a massless parton with momentum p l splits into massless daughter partons with momentap l andp m+1 , we define (Ref. [4], Eq. (4.19)) where Q is the total momentum of the final state partons before the splitting. For an initial state splitting in which a massless initial state parton with momentum p l (l = a or b) splits in backward evolution into a new massless initial state parton with momentump l and a massless final state parton with momentump m+1 , we define (Ref. [10], Eq. (4.1)) In addition to y, the splitting functions in H I (t) depend on a momentum fraction z (Ref. [13], Eqs. (7) and (9)). Various definitions of z are possible. For a final state splitting, Deductor useŝ where the auxiliary lightlike vectorñ l is defined using the total momentum Q of the final state partons: For the splitting of initial state parton "a" from hadron A, z is the ratio of momentum fractions η a before andη a after the splitting: Here p b is the momentum of the initial state parton from hadron B. There is a third splitting variable, the azimuthal angle φ ofp m+1 about the direction of the mother parton, p l in the rest frame of Q. We will denote the splitting variables {y, z, φ} collectively by ζ p and the flavor choice in the splitting by ζ f . We denote {ζ p , ζ f } by ζ, as described in more detail in Appendix A.
The default choice of the shower ordering variable in Deductor is Λ 2 , defined by where Q 0 is the total momentum of the final state particles at the start of the shower. (See Ref. [13], Eq. (5)). Then the shower time t is the logarithmic variable (Ref. [13], Eq. (A.6)), so that y → 0 implies t → ∞.
The splitting functions in H I (t) are singular in the limit y → 0. There are two kinds of singularities. There are collinear singularities, corresponding to y → 0 at fixed z. There are also fixed angle, soft singularities, corresponding to emission of a soft gluon carrying momentum fraction (1 − z): y → 0 and z → 1 with fixed y/(1 − z). There are also collinear×soft singularities, in which y → 0, z → 1, and y/(1 − z) → 0. We recall [7] that the LC+ approximation is exact for collinear splittings and for collinear×soft splittings. It it approximate only for fixed angle soft splittings.
In this section we explore the contribution of the collinear and collinear×soft regions compared to the contribution of the fixed angle soft contribution to the operators ∆H(τ ) and X (τ, t 0 ). We do this by decomposing the operators involved into parts that get contributions only from the fixed angle soft region and parts that get contributions from everywhere else. The Deductor code is not organized using this decomposition, but nevertheless the analysis of this section can help us to understand the behavior of the approximations that we use to go beyond the LC+ approximation.
We can divide H I (t) into two contributions [4], The contribution H coll. (t) contains the collinear and also the collinear×soft singularities, while H soft (t) contains only the fixed angle, soft singularities. We review the definitions for this decomposition in Appendix. B. Here we state only a few important properties of the two contributions.
The contribution H coll. (t) comes from the square of a Feynman graph in a physical gauge in which the new parton m + 1 is emitted from a given parton (e.g. l or a), as illustrated in Fig. 1. This gives a function whose y → 0 limit is a DGLAP splitting function P (z), including its 1/(1 − z) singularity.
The contribution H soft (t) comes from interference of two emissions: a soft gluon emitted from one parton in the ket state and the same soft gluon emitted from a different parton in the bra state, as illustrated in Fig. 2. As we will see in Appendix B, in a physical gauge, this contribution has at most integrable singularities in the limit in which the emitted gluon becomes collinear to either of the two emitting partons [5]. Soft real emission from parton l interfering with emission from parton k.
In Eq. (18), we defined another decomposition of H I (t), According to the definition of the LC+ approximation in Ref. [7], the LC+ approximation is exact for splittings corresponding to the square of a Feynman graph in a physical gauge. We can divide H coll.

I
(t) into its LC+ approximation and a remainder, as in Appendix B, We have That is, ∆H(t) has only a fixed angle, soft contribution, but has no collinear contribution. On the other hand, the LC+ approximation is not exact for H soft (t), so that in the decomposition both contributions are non-zero. The Sudakov exponent V(t) can be decomposed in an analogous manner, Here V iπ (t) is the contribution from the imaginary part of virtual graphs and is given in Eq. (9). The splitting functions V coll. (t) and V soft (t) represent splittings that did not happen. They are the same as those in H I (t) except for their color structure and except that we now integrate over z and φ. Thus if we write we have In this paper, we do not introduce an LC+ approximation for V iπ (t), so as in Eq. (19). The corresponding color structures [7] are illustrated in Fig. 3 for V coll. (t) and Fig. 4 for V soft (t). The color structure of V soft (t) is non-trivial, so that we cannot exponentiate V soft (t) in any simple way. Fortunately, the color structure of V coll. (t) is very simple: for gluon emission from a quark, it is proportional to a unit operator on the color space times a factor C F , for gluon emission from a gluon, it is proportional to C A , and for a g → q +q splitting, it is proportional to T R = 1/2 .
The simplicity of the color structure of V coll. (t) has an important consequence: even if the color state to which we apply V coll. (t) changes, applying V coll. (t) to this state still returns the state times C F , C A , or T R depending on the flavor of the splitting parton. This means that V lc+ coll. (t) commutes with V lc+ soft (t) and the full ∆V(t). Because V lc+ coll. (t) commutes with V lc+ soft (t), N lc+ (t 2 , t 1 ) in Eq. (16) takes the form Additionally, in Eq. (26) for X (t, t 0 ), all of the factors of exponentials of V lc+ coll. (τ ) commute with the other factors in X (t, t 0 ), giving That is, only the soft contribution to V lc+ (τ ) contributes to X (t, t 0 ). Furthermore, according to Eq. (43), only V iπ (τ ) and the soft contribution to V(τ ) contributes to ∆V(τ ). Thus everywhere in Eq.
We note that in calculating X (t, t 0 ) according to Eq. (47), we exponentiate V lc+ soft (τ ) and expand perturbatively in powers of ∆V soft (τ ) and V iπ (τ ). We also note that the operator V coll. (τ ) carries both collinear and collinear×soft singularities and thus can contribute double large logarithms for each power of α s for some observables. Since this operator does not appear at all in X (t, t 0 ), the factor X (t, t 0 ) in Eq. (27) gives only single logarithms for each power of α s .

VII. COLOR STATES
In the previous versions of Deductor, we always evolved a single basis state and after every step of the shower we always expanded the color state on basis states. The corresponding sums over basis states were then performed by Monte Carlo summation: picking one term at random. This is useful in the implementation since it keeps the code rather simple. However when we consider the effect of the fixed angle soft radiation, we can have serious numerical problems because performing all sums over color basis states by Monte Carlo summation leads to greater fluctuations than one wants. Instead, we can always select a unique basis state in the momentum and flavor space but in the color space we can use a linear combination of the color basis states. For this purpose, we define a statistical state with definite momentum and flavor choice but a more general color state as The state is labeled by the function ψ giving the coefficients of the expansion of the color state in color basis states.

VIII. THE LC+ NO SPLITTING OPERATOR
From Eq. (27), we see that the no-splitting operator in the LC+ approximation plays an important role in shower evolution. The LC+ approximation is defined in such a way that every basis state {p, f, c , c} m is an eigenstate of N lc+ (t 2 , t 1 ) [7], When the LC+ no-splitting operator acts on a generic state we have The result is a linear combination of the basis vectors, as in {p, f, ψ} m , but the terms are weighted by the corresponding Sudakov factor. In this section, we rewrite the LC+ no-splitting operator as the product of an average Sudakov factor ∆(t 2 , t 1 , {p, f, ψ} m ) for the generic color state ψ and a weight factor given by an operator Φ. We can use the average Sudakov factor to select splitting variables.
(52) 6 The function The Sudakov exponent has the form 7 Recall from Sec. V that ζ = {ζ p , ζ f } is a shorthand notation for the splitting variables in momentum and flavor. Then dζ stands for integrating over the momentum splitting variables and summing over the new flavors as in Eqs. (A16) and (A17) below. The function T specifies the shower time according to eq. (35). The function λ lk ({p, f } m , ζ) is a rather complicated non-negative color-independent function given in Eq. (A4). The function χ(k, l, {c , c} m ) is a simple momentum-independent function We want to define an averaged Sudakov exponent that is simple in color and has all the soft-collinear and collinear singularities. This can be done many ways but in this paper we define this via an averaged characteristic function, ξ(k, l, {ψ} m ). In general it can depend on the {ψ} m state. In order to recover the exact collinear and soft-collinear singularities in the averaged Sudakov exponent we should have We also require We give two possible definitions of the averaged characteristic function ξ(k, l, {ψ} m ) below. Using whatever definition is chosen, we can define an averaged Sudakov exponent 7 The function λ lk ({p, f }m, ζ) χ(k, l, {c , c}m) was called λ({p, f, c , c}m, l, k, {p,f } m+1 , )) in Eq. (7.6) of Ref. [7], although the c argument seems to be missing in Ref. [7]. and the corresponding Sudakov factor Now the no-splitting operator is the product of the new Sudakov factor ∆(t 2 , t 1 , {p, f, ψ} m ) and an operator Φ that supplies a correction factor: where the operator Φ(t 2 , t 1 ; ψ) is defined as This definition gives us Our expectation is that the difference between λ({p, f, c , c} m , τ ) and λ({p, f, ψ} m , τ ) will be reasonably small, so that the weight factor created by Φ(t 2 , t 1 ; ψ) will be close to 1. The reason for this expectation is that splitting functions in λ({p, f, c , c} m , τ ) are independent of the color state {c , c} m in the limit of collinear and collinear×soft splittings. This means that the difference between λ({p, f, c , c} m , τ ) and λ({p, f, ψ} m , τ ) in Eq. (62) is sensitive only to fixed-angle soft splittings, which we expect are not numerically very important. In designing Deductor 3.0.0, we considered and implemented two version of the averaged characteristic function. Define a function p(ψ, {c , c} m ) that assigns probabilities to the basis states |{c , c} m ). The probabilities are positive and normalized to The choice of the probabilities is largely arbitrary and the physical quantities are independent of them. One choice tries to emphasize the importance of the basis state, (64) Here I 0 ({c , c} m ) is defined using the color overlap of color basis states using the U (N c ) group instead of SU (N c ): This choice of the probability function emphasizes the color basis states in ψ that don't have big I 0 . Now, with these probabilities, the first averaged characteristic function is defined as The other choice is much simpler. It doesn't depend on the color state at all but depends only on the flavor and number of the partons. We define (67) Both choices of the characteristic function have some advantages and disadvantages. We will discuss them in the next section.

IX. EVOLUTION
We are now in a position to rewrite the evolution equation Eq. (27) in a computationally useful form. We start by simply applying Eq. (27) to a state {p, f, ψ} m : Next, we use Eq. (60) for N lc+ and, inside the integral, we multiply and divide by the normalization factor λ({p, f, ψ} m , τ ). This gives us In order to use this, we need to define the color coefficients of X Φ ψ : Then Now we examine what the splitting operator H I (τ ) = H lc+ (τ ) + ∆H(τ ) does when applied to an arbitrary basis state, When k = l, only the LC+ term is present: Now we factor k =l λ lk ({p, f } m , ζ) ξ(k, l, {ψ} m ) from the second term. Then we have Inserting this into Eq. (71) and using ξ(l, l, {ψ} m ) = 2, we obtain The color matrix X transforms the original color state ψ for m partons to a linear combination of color basis states {c , c} m under the color transformation provided by the operators X and Φ. Then the expression under the square brackets is the color matrix that transforms the color state {c , c} m to {ĉ ,ĉ} m+1 . This matrix depends on the hard state {p, f } m , the splitting variables ζ, and the parton labels l, k.
We now have a form that is useful for calculations. The factors λ∆ dτ in the integral over τ say that we should pick a next shower time τ with a probability determined by the total color state ψ. This is different from the probability that we would have with a single color basis state {c , c} m . However, we account for the probabilities for each basis state in the factor Having fixed τ , we also need to pick l, k, and the split-ting variables ζ. We pick l, k, and ζ with probability According to Eq. (58), this gives us a properly normalized probability, lk dP l,k,ζ = 1.
We are left with a sum over colors {ĉ ,ĉ} m+1 for the m + 1 partons after the splitting. To perform this sum, we must chart a sensible course [22]. At one extreme is Scylla: we could perform this sum Monte Carlo style, picking one color state {ĉ ,ĉ} m+1 at random and accumulating a weight factor given by the remaining factors in Eq. (76). This leads to large fluctuations that render the calculation useless. At the other extreme waits Charybdis: we could accumulate all of the color states {ĉ ,ĉ} m+1 into a new combined color state ψ. This leads to color states ψ that contain so many basis states {ĉ ,ĉ} m+1 that the calculation crashes. We choose a middle course. The algorithm can then be adjusted to optimize performance. Its details are not important, at least conceptually.
The algorithm described above is implemented in Deductor with both ξ 1 (k, l, {ψ} m ) and ξ 2 (k, l, {ψ} m ) averaged characteristic functions. We also implemented an alternative way to do the shower evolution. It also uses the ξ 1 (k, l, {ψ} m ) characteristic function and the difference is only in the treatment of the H I (t) operator. In the alternative shower algorithm for the splitting operator we have This way we have a little bit faster algorithm but we have larger fluctuations from the second term under the square bracket due to the potentially small denominator.

X. APPROXIMATIONS
We now need some approximations in order to have an algorithm for color evolution that can be executed on a finite size computer in a finite amount of time. The key requirement is that one should be able to control the level of approximation, so that one can obtain more nearly exact results if one has greater computer resources available.
To start, we put a limit on the color suppression index I in the shower evolution operator U(t, t 0 ). As defined in Ref. [7], the color suppression index is obtained from two factors. The first is the number p E of explicit powers of 1/N c that arise from choosing the color suppressed term in the Fierz identity for g → q +q splittings that have occurred in the shower history so far. The second is the number of powers of 1/N c in the overlap of the bra and ket color states using U(N c ) color: (80) The cross section at the end of the shower will be suppressed by a color factor that is at least as small as 1/N I c [7], so there is little point in keeping contributions with large I. Thus we use the full SU(N c ) evolution for color states for which I − I hard ≤ I max , where I hard is the color suppression index of the hard scattering state at the start of the shower and I max is a parameter that we choose.
A sensible choice for I max is I max = 4. We can improve the approximation that results from limiting I by increasing I max , but there is a cost. Increasing I max increases statistical fluctuations in the results after a fixed amount of computer time, so that we need more computer time to achieve the same statistical accuracy.
If the shower operator reaches I − I hard ≥ I max , it switches to an approximate shower based on the color group U(N c ) instead of SU(N c ). We also omit any further contributions from ∆H and ∆V. Thus contributions proportional to 1/N I c are calculated only approximately. Next, we put a limit on the number of times, N thr. ∆ , that the operator ∆V = ∆V Re + V iπ defined in Eqs. (18) and (19) appears in the threshold operator U V . This operator gives results as an expansion in powers of ∆V and we retain only those terms with no more than N thr. ∆ factors of ∆V.
Then, we put limits on the number of times that the operators ∆H, ∆V Re and V iπ appear in the shower evolution operator U(t 2 , t 1 ). As described in Sec. IX, U(t 2 , t 1 ) produces terms proportional to [ We choose parameters N ∆ , N Re and N iπ . We retain only terms with A + B ≤ N Re and C ≤ N iπ and gives more accurate results in the limit of long computer running times and large computer memory. However, it could use more memory than is available and it increases statistical fluctuations in the results after a fixed amount of computer time, so that we need more computer time to achieve the same statistical accuracy.
We can prescribe different accuracy parameters N Re , N iπ , N ∆ and I max to successive shower time intervals. One would do this if one suspected that, for the observable being measured, the first splitting steps of the shower are more important than later steps, with softer splittings. In the simplest application, we assign N Re , N iπ , N ∆ , and I max to the interval t (0) < t < t (1) . Here t (0) is the starting time of the shower as determined by the hard interaction that initiates the shower and t (1) > t (0) , perhaps t (1) = t (0) + 5. Then we can use an LC+ shower (N ∆ = 0) with the same I max until the shower ends. The shower ends because the splitting kernel has a cutoff built into it that stops splittings at a lower cutoff for the transverse momentum in a splitting. We choose this cutoff to be k min T = 5 GeV.
The Deductor 3.0.0 user interface also allows one to specify accuracy parameters for successive shower intervals determined by a fixed number of splitting steps.

XI. PROBABILITY CONSERVATION
Version 3.0.0 of Deductor implements the algorithm outlined above. It is capable of producing cross sections, as we will see below in Sec. XII. However, it is not so easy to check whether it is producing correct cross sections including non-leading color effects, as intended.
In this section, we present a test of the inner workings of the program by testing the probability conservation property (8), Here ρ(t 0 ) is an arbitrary statistical state at shower time t 0 . Then 1 ρ(t 0 ) is the total cross section measured for that state, obtained by integrating the differential cross section over all parton variables and taking the trace of the color density matrix. After generating a shower with U(t, t 0 ), we have a much more complicated statistical state with typically many more partons. We then measure the total cross section associated with this state. The total cross section should be the same.
Note that it is a simple consequence of the shower evolution equation (17) that probability conservation, Eq. (81), holds in the LC+ approximation and also with the inclusion of contributions from the operators ∆H, ∆V Re , and V iπ . The evolution equation (17) leads us to the evolution equation (27). In this form of the evolution equation, it is no longer self-evident that probability conservation works, since ∆H and ∆V are treated very differently now. With further manipulations, we arrive at the algorithm implemented in Deductor, Eq. (76). Now, probability conservation must still work but it is not a property that one would guess from examining Eq. (76) without knowledge of its derivation. Thus it can be a powerful test of Eq. (76) and its implementation in code to see whether probability conservation (81) works within the statistical accuracy of the calculation.
We will check whether for each value of t 0 . We will expand the left-hand side of (82) in powers of ∆H, ∆V Re , and V iπ and examine the terms proportional to The term with A = B = C = 0 gives us the LC+ approximation to U(t, t 0 ), which obeys Eq. (82). Thus the term with A = B = C = 0 gives us the 1 on the right hand side of Eq. (82). Consequently, the other times must combine to give zero. Furthermore, if we were to replace ∆H −∆V Re by λ Re (∆H − ∆V Re ) and V iπ by λ iπ V iπ , the relations 1 (∆H − ∆V Re ) = 0 and 1 V iπ = 0 tell us that Eq. (82) holds order by order in λ Re and λ iπ . Thus all of the contributions with fixed values of A + B and C, other than A = B = C = 0, must sum to zero.
In this test, the shower is limited to the shower time interval t 0 < t < t 0 + 5. We choose the maximum color suppression index to be I max = 4. We also limit the number of ∆H and ∆V operators as specified below.
In our test, ρ(t 0 ) is the state produced by a 2 → 2 hard scattering. The two final state partons have absolute value of transverse momenta p 1,T = p 2,T = p T and c.m. energy squared (p 1 + p 2 ) 2 =ŝ. With our shower time definition and the choice of the starting time for the shower from Ref. [13], the starting shower time is t 0 = log(4ŝ/(9p 2 T )). Thus t 0 = log((8/9)(1 + cosh(y 1 − y 2 )). We generate hard scatterings with a wide range of transverse momenta, from P T = 20 GeV to P T = 5 TeV. The smallest possible value for t 0 is about 0.6 and, with small values of P T , t 0 can range up to more than 10. We generate a range of t 0 values and plot results in bins of t 0 , so that t 0 has a definite value in Eq. (82).
By default, Deductor has an operator U V that comes between the hard state ρ(t 0 ) and the probability preserving shower operator U(t, t 0 ). This operator sums threshold logarithms and thereby changes the cross section. For this test, we turn U V off.
We wish to test Eq. (82) in the presence of approximations with respect to color. However, we caution that exact agreement with Eq. (82) is not to be expected, since there are small systematic errors within the Deductor calculation that come from sources other than limits on the color treatment. For instance, there are inevitably approximations in our use of the parton distribution functions. Additionally, there are some functions defined as integrals that are too complicated to be performed analytically; for these, we use gaussian numerical integration. We believe that these systematic errors are smaller than 1%, but we have not systematically checked their size.
We will test Eq. (82) in three ways. First, we turn off V iπ and use only ∆H and ∆V Re . Then we turn off ∆H and ∆V Re and use only V iπ . Finally, we use ∆H and ∆V Re and V iπ together.
We now try the first test in which we turn off V iπ and consider up to four insertions of ∆H and ∆V Re . That is, we set N Re = 4, N iπ = 0, and N ∆ = 4.
In With A = B = 0, we are looking at probability conservation from the LC+ approximation. We see that this contribution is 1 to within small statistical fluctuations.
Here and in later graphs, we plot error bands that represent the estimated statistical error in the Monte Carlo data. We can also assess the statistical errors by looking at the fluctuations from one bin in the plot to the next. In Fig. 5, the statistical errors are hardly visible, but they are more visible in later plots.
In Fig. 6, we look at the contribution to the left hand side of Eq. (82) from contributions with A + B = 1. We plot separately the contributions from ∆H and ∆V Re along with their total. We see that these contributions are typically between ±1% and ±4%. However, they cancel, giving a total that is smaller than about ±0.2%. It is remarkable to us that these contributions cancel to the extent seen in the figure since the methods of calculation for ∆H and ∆V Re are very different. In Fig. 7, we look at the contribution to the left hand side of Eq. (82) from contributions with A + B = 2. We plot separately the contributions from ∆H 2 , ∆H ∆V Re , and ∆V 2 Re , along with their total. We see that the individual contributions are typically ±2%. However the total is smaller than about ±0.2%.
In Fig. 8, we look at the contributions with A + B = 3. We plot separately the contributions from ∆H 3 , ∆H 2 ∆V Re , ∆H ∆V 2 Re , and ∆V 3 Re , along with their total. The individual contributions are small, of order ±10 −3 . This small size is understandable because each contribution is proportional to α 3 s . The statistical fluctuations are almost as large as the individual contributions. The total is smaller than about ±10 −3 , certainly smaller than the sum of the absolute values of the individual contributions.
Finally, in Fig. 9, we look at the contributions with A + B = 4. We plot separately the contributions from ∆H 4 , ∆H 3 ∆V Re , ∆H 2 ∆V 2 Re , ∆H ∆V 3 Re , and ∆V 4 Re , along with their total. The individual contributions are again small, of order ±5 × 10 −4 . The contributions with the most factors of ∆H have fluctuations roughly as large as the contributions. The total inherits the fluctuations of the individual contributions but is consistent with zero within its fluctuations.
We see that all contributions proportional to [∆H] A [∆V Re ] B with A + B = 1, 2, 3, or 4 cancel to within their statistical accuracy, confirming Eq. (82) for these contributions.
We now try the second test of Eq. (82), in which we turn off ∆H and ∆V Re and consider up to four insertions of V iπ . That is, we set N Re = 0, N iπ = 4, and N ∆ = 4.
We divide the left hand side of Eq. (82) into parts that we can plot separately. The operator V iπ in Eq. (9) is the sum of two terms, V iπ = V L iπ + V R iπ where in V L iπ the color operator T a · T b operates on ket color states and in V R iπ the color operator T a · T b operates on bra color states. Acting on a color basis state {ĉ} m , T a · T b produces terms proportional to new basis states {ĉ} m . When partons "a" and "b" are color connected in {c} m , one of these new basis states is a constant times the original basis state {c} m . Thus we can write V L iπ = V L iπ + I L iπ , where I L iπ is the part of V L iπ that returns a constant times the original color basis state and V L iπ is the part of V L iπ that changes the color basis state. We apply the same decomposition to V R iπ . As a bookkeeping measure, it is convenient to decompose V iπ as where I iπ = I L iπ + I R iπ . That is, we label contributions according to whether the ket color state was changed, V L iπ , or the bra color state was changed, V R iπ , or the color state was unchanged, I iπ . In the figure legends, we abbreviate these operators as L, R, and I.
In Fig. 10, we look at the contribution to the imaginary part of the left hand side of Eq. (82) from contributions proportional to V iπ . We plot separately the contributions from V L iπ , V R iπ , and I iπ along with their total. We see that the individual contributions are of order ±4% and that these contributions cancel, giving a total that is smaller than about ±0.2%.
In Fig. 11, we look at the contribution to the real part of the left hand side of Eq. (82) from contributions proportional to V 2 iπ . We plot separately the contributions along with their total. We see that some of the individual contributions are of order ±1. It is to be expected that the contributions from powers of V iπ are larger than those from the same powers of ∆H and ∆V Re because α s in V iπ comes with a factor 4π. The contributions from these terms in V 2 iπ cancel, giving a total that is smaller than about ±0.2%.
In Fig. 12, we look at the contribution to the imaginary part of the left hand side of Eq. (82) from contributions proportional to V 3 iπ . We plot separately the con- , along with their total. We see that some of the individual contributions are of order ±0.3. Again, the contributions from these terms in V 3 iπ cancel, giving a total that is smaller than about ±0.2%.
In Fig. 13, we look at the contribution to the real part of the left hand side of Eq. (82) from contributions proportional to V 4 iπ . We plot separately the contributions , along with their total. We see that some of the individual contributions are of order ±0.4. Yet again, the contributions from these terms in V 4 iπ cancel, giving a total that is smaller than about ±0.2%.
We see that all contributions proportional to (V iπ ) C with C = 1, 2, 3, or 4 cancel to within their statistical accuracy, confirming Eq. (82) for these contributions.
Finally, we consider up to four insertions of ∆H, ∆V Re When we examine contributions to the real part of the cross section with A + B + C = 1, the contribution from V iπ is not present because this contribution in imaginary. Thus we obtain the graph shown in Fig. 6, in which the needed cancellations work up to statistical fluctuations, which are smaller than about ±0.2%.
In Fig. 14, we look at contributions with A + B + C = 2. We plot separately the contributions from ∆H 2 , ∆H ∆V Re , ∆V 2 Re , and V 2 iπ along with their total. We do not break the V 2 iπ contribution into separate parts as we did in Fig. 11. Recall from Fig. 11 that the separate parts of the V 2 iπ contribution are of order 1. We see that the total with everything included vanishes up to statistical fluctuations, which are smaller than about ±1% and are dominated by the statistical fluctuations in V 2 iπ . In Fig. 15, we look at contributions with A + B + C = 3. We plot separately the contributions from ∆H 3 , ∆H 2 ∆V Re , ∆H ∆V 2 Re , ∆V 3 Re , ∆HV 2 iπ and ∆V Re V 2 iπ along with their total. We see that the total vanishes up to statistical fluctuations, which are smaller than about ±1% and are dominated by the statistical fluctuations in the terms that contain V 2 iπ . Finally, in Fig. 16, we look at contributions with A + B + C = 4. We plot separately the contributions from ∆H 4 , ∆H 3 ∆V Re , ∆H 2 ∆V 2 Re , ∆H ∆V 3 Re , ∆V 4 Re , ∆H 2 V 2 iπ , ∆H ∆V Re V 2 iπ , ∆V 2 Re V 2 iπ , and V 4 iπ , along with their total. In this plot, the contributions with no powers of V iπ are small. We see that the total vanishes up to statistical fluctuations, which are smaller than about ±1% and are dominated by the statistical fluctuations in the terms that contain powers of V iπ . In summary, we have seen in some detail that probability conservation, Eq. (82), works.

XII. A SAMPLE CROSS SECTION
Version 3.0.0 of Deductor, as described above, can be used to calculate cross sections with color treated beyond the LC+ approximation. We leave an examination of the phenomenology of color corrections to a future work. However, in this section we demonstrate at least that Deductor with non-leading color effects can be used to calculate a physical cross section. We choose a cross section for which we expect that color effects beyond the LC+ approximation will be small, namely the one jet inclusive cross section dσ/dP T for jets with rapidities in the range −2 < y J < 2 as a function of the jet transverse momentum P T for proton-proton collisions at √ s = 13 TeV. For the renormalization and factorization scales in the Born cross section with which the shower begins, we choose µ R = µ F = P Born T / √ 2. The shower cross section depends on the Λ-ordering starting scale for the shower, which we choose to be µ s = (3/2)P Born T , which corresponds to a shower time t 0 = log(4ŝ/(9[P Born T ] 2 )). We have examined this cross section in some detail in Ref. [13] in the LC+ approximation. Here we limit our t0 Im σ(showered)/σ(Born) investigation to the effects of color beyond the LC+ approximation.
The shower stops when the transverse momentum in a splitting reaches a cutoff k min T . In this section, we choose k min T = 5 GeV. In Ref. [13], we used k min T = 1 GeV for final state splittings and k min T = 1.295 GeV for initial state splittings. The 5 GeV choice has the advantage of making the code run faster. If we wanted to correct to the choice in Ref. [13], we would apply a correction factor that ranges from about 0.93 at P T = 300 GeV to 1.0 at P T = 3.5 TeV. We could also apply a correction factor of about 0.98 to account for non-perturbative effects [13].
The formula that represents what Deductor does is the following, taken from Eq. (134) of Ref. [14] with some simplification of the notation: The statistical state ρ hard includes the scattering matrix elements and a factor containing the parton distribution functions for the two incoming partons for a range of transverse momenta and rapidities in the hard scattering. In principle, ρ hard should include next-to-leading order (NLO) corrections with their accompanying subtractions, as described in Ref. [14]. However, we have not implemented the NLO corrections in Deductor 3.0.0. This limits the accuracy in the calculation. The initial shower time t 0 is determined from the kinematics of the initial hard scattering: t 0 = log(4ŝ/(9[P Born T ] 2 )). Thus t 0 in Eq. (84) is really an operator that gives log(4ŝ/(9[P Born T ] 2 )) as its eigenvalue for a statistical basis state in the expansion of ρ hard .  The operator U V (t f , t 0 ) implements the approximate summation of threshold logarithms, as described in Refs. [13] and [14]. We have, however, dropped some numerically unimportant terms from U V (t f , t 0 ) compared to Ref. [13]. The next operator, U(t f , t 0 ), generates the shower. The shower stops when the transverse momentum in a splitting is smaller than a cutoff value k min T = 5 GeV. Then the operator O J specifies the jet measurement that we want to make. Finally, the bra state 1 is an instruction to integrate and sum over all of the parton variables, including taking the trace over the color variables.
Since dσ/dP T is a steeply falling cross section, we display results for the ratio to the next-to-leading order cross section [23], In the plots that follow, we choose the maximum color suppression index to be I max = 4. This applies both within the LC+ approximation and when operators ∆H and ∆V are allowed.
In the upper red curve in Fig. 17, we plot K(P T ) with two units of color beyond the LC+ approximation in the shower and one unit of extra color in the threshold factor: N ∆ = N Re = N iπ = 2 and N thr. ∆ = 1. In the shower, the inclusion of extra color applies in the first 5 units of shower time, t 0 < t < t 0 + 5. After that, we use an LC+ shower. The upper blue curve is K(P T ) in the LC+ approximation, N ∆ = N Re = N iπ = 0 and N thr. ∆ = 0. We label these as "Full" calculations, indicating that they include the threshold factor U V . The threshold factor is quite important: in the lower, dashed red and blue curves, we plot K(P T ) with U V turned off, so that we have only the standard ("Std."), probability preserving shower generated by U(t f , t 0 ). The dashed red curve is calculated with N ∆ = N Re = N iπ = 2 in the first 5 units of shower time while the dashed blue curve is calculated in the LC+ approximation, N ∆ = N Re = N iπ = 0.
We see that there is a small effect from going beyond the LC+ approximation both with and without the threshold factor. This approximately 3% effect appears both in the Full calculation and the Std. calculation.
We find that the small effect from extra color comes from the V iπ operators. To see this, we plot in Fig. 18 the same curves as in Fig. 17 but without V iπ . Now there is no difference, within statistical fluctuations, between the red and blue curves.
It is understandable that V iπ affects the jet cross section. We have seen that V iπ , with its factor of 4πα s , is not a small operator. Thus contributions from V 2 iπ in Fig. 11 are of order 0.2. These contributions cancel when we use a completely inclusive measurement operator, as we saw in Fig. 11. However the different contributions have different color states, so that we can radiate more or fewer gluons out of the jet cone when we use a jet measurement operator. This can change the jet cross section.
What happens if we add more powers of ∆H and ∆V? We can double the amount of extra color by using N thr. and ∆V in U(t f , t 0 ) changes the cross section by a factor of only about 1.03. Thus we might expect that adding two more powers of ∆H and ∆V will further change the cross section by a factor of only 1 + (0.03) 2 ≈ 1.001. However, this expectation could be wrong. Perhaps an unanticipated effect will change the cross section by a much larger factor, say 1.1. We can check by simply doing the calculation. Fig. 19, we plot results as in Fig. 17 but with N ∆ = N Re = N iπ = 4 and N thr. ∆ = 2. The red curves show results with the extra color insertions while the blue curves show the results with just the LC+ approximation. The upper, solid curves are for the Full calculation including the threshold factor while the lower, dashed curves are for the Std. calculation with only the probability preserving parton shower. We see that the statistical fluctuations are much larger than they were with just two units of extra color. 8 We now cannot see the approximately 3% change in the cross section that resulted from adding two units of extra color. However, if four units of extra color added 10% to the cross section, we could see that even with the greater statistical fluctuations.

XIII. CONCLUSIONS
We can offer three arguments for pursuing corrections to a leading order parton shower beyond the leading color approximation.
When one uses an NLO calculation for the hard scattering that initiates a shower, the hard scattering calculation should be matched to the parton shower. (Deductor does not currently include this matching.) The NLO calculation naturally includes contributions beyond the leading order in 1/N 2 c . It is not impossible to perform the matching when the shower lacks these contributions. However, it is clearly best if the shower includes the same color accuracy as the NLO hard scattering calculation for the color density matrix.
One can also hope to eventually have a parton shower algorithm in which the splitting kernels are correct to order α 2 s instead of just α s . However, 1/N 2 c is of roughly the same size as α s . Thus it would seem as important to include 1/N 2 c corrections to the order α s splitting functions as it is to include α 2 s contributions to the shower splitting functions.
Perhaps most importantly, there may be processes in which corrections beyond the leading order in 1/N 2 c are numerically important because these corrections multiply large logarithms [16,[24][25][26][27][28][29][30][31][32][33][34]. A parton shower is a promising way to investigate these effects, but evidently the shower must incorporate (1/N 2 c ) k corrections. Now, we briefly review the argument of this paper. We view a parton shower as an application of the renormalization group, proceeding from harder interactions to softer interactions. We take the proper framework for the shower to be quantum statistical mechanics. Simple classical statistical mechanics is not sufficient for two reasons. First, we need to account for quantum interference between emission of a soft gluon from one parton and emission from another gluon. Second, parton color is a quantum degree of freedom. Using quantum statistical mechanics for color requires one to consider the {c} m {c } m density matrix in color. (The density matrix in parton spin space is required also, but in this paper we ignore spin.) For a lowest order parton shower, we need two operators, which we call H I (t) and V(t). H I (t) describes parton splittings, which increase the number of partons by one. It has a certain color structure, which one simply reads off from the Feynman rules for QCD. V(t) leaves the number of partons unchanged. It represents an approximation to one loop virtual graphs. It has a certain color structure, which one simply reads off from the Feynman rules for QCD.
The operators H I (t) and V(t) are related by the equation, in the notation of this paper, 1 H I (t) = 1 V(t). This means that the ideal shower based on H I (t) and V(t) preserves probabilities. Fundamentally, this property arises from the fact the infrared divergences of QCD cancel between real and virtual graphs, which in turn follows from the fact that the quantum evolution operator U (t 2 , t 1 ) is unitary [35].
We know the color structure of H I (t), so it is straightforward to incorporate H I (t) into computer code, as in Refs. [15,17]. However, we need also to incorporate V(t) in order to be consistent with quantum mechanics.
It is, unfortunately, not easy to incorporate V(t) into computer code for a parton shower. The reason is that, in the traditional formulation of a parton shower, one needs a Sudakov factor consisting of the exponential of an integral of V(t), but this exponential is an operator on the color space. As the number of partons increases, the dimensionality of the color space becomes enormous, so the exponential of an integral of V(t) becomes difficult to calculate.
One commonly applies the leading color (LC) approximation in a practical parton shower. Here H I (t) is approximated by an operator H lc (t) and V(t) is approximated by an operator V lc (t). This approximation gets cross sections right to leading power in 1/N 2 c , where N c = 3 is the number of colors. It has the property that V lc (t) is diagonal in color, so that is easy to exponentiate.
The LC approximation also has the crucial property 1 H lc (t) = 1 V lc (t). For this reason the LC shower preserves probabilities.
Deductor uses an improved approximation, the LC+ approximation, as its starting point for treating color. This approximation retains some contributions that the LC approximation drops. In this approximation H I (t) is approximated by an operator H lc+ (t) and V(t) is approximated by an operator V lc+ (t). The approximation has the property that V lc+ (t) is diagonal in color, so that is easy to exponentiate. It also has the property that 1 H lc+ (t) = 1 V lc+ (t), so that the LC+ shower pre- serves probabilities.
The LC+ approximation becomes exact in the limit of collinear emissions or soft×collinear emissions. That is, the only singular limit for emissions in which the LC+ approximation is not exact is the limit of fixed angle soft emissions. This feature, which is not shared by the LC approximation, is important for the working of the algorithm described in this paper.
Deductor is organized as a dipole shower, in which there is quantum interference between emission of a gluon from a parton with label l and another parton with label k. The symmetry between l and k is removed in such a way that there is a singularity when the emitted gluon becomes collinear with parton l but not with parton k. The momentum mapping used in Deductor depends on the choice of l but not the choice of k. This feature is important in the algorithm used in this paper.
We can retain terms with up to N ∆ powers of ∆H and ∆V in the first t ∆ units of shower time and we can retain N thr.
∆ powers of ∆V in the threshold operator U V . We also impose a limit I max on how large the color suppression index I, Eq. (80), can grow. This limits the number of powers of 1/N c that we keep.
The resulting algorithm for a parton shower is approximate with respect to color but the approximation is systematically improvable by making N ∆ , N thr. ∆ , t ∆ , and I max larger, at the cost of requiring more computer memory or more computer time to reach the same level of Monte Carlo statistical accuracy. This algorithm is implemented in public computer code. 9 We have tested whether the cancellations between ∆H and ∆V actually work so as to produce a probability preserving shower. Within the accuracy of the calculations, these cancellations do work.
We have not used the new version of Deductor to investigate cross sections in which color beyond the LC+ approximation might play an important role. We expect to carry out such investigations in future work. However, we have calculated the one jet inclusive cross section beyond the LC+ approximation just to check how well the program works in calculating a physical cross section.
We now turn to the outlook for future work. We expect that the algorithm presented here will not be the last word in algorithms for this purpose. Surely it is possible to do better. Indeed,Ángeles Martínez, De Angelis, Forshaw, Plätzer, and Seymour [16] have provided a formalism for the description of soft gluon emissions that is similar in some ways to the general formalism [4,7] on which this paper is based. If the approach of Ref. [16] can be extended to include the collinear singularities of QCD, then it will be of great interest to see if there can be a practical implementation of the resulting formalism. Perhaps such an implementation will be able to outperform what this paper provides.
Our numerical investigations suggest that V iπ is effectively a larger operator than V Re . For this reason, it may be better to include V iπ in V lc+ instead of ∆V. This means that one would need to numerically exponentiate V iπ . This will cost computer resources, so it remains to be seen if this is a better option.

ACKNOWLEDGMENTS
This work was supported in part by the United States Department of Energy under grant DE-SC0011640. This project has benefited from the participation of the authors at the Munich Institute for Astro-and Particle Physics program "Automated, resummed and effective: precision computations for the LHC and beyond" as well as 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. DS thanks the Erwin Schrödinger Institute of the University of Vienna for its hospitality at its program "Challenges and Concepts for Field Theory and Applications in the Era of the LHC Run-2." This work benefited from access to the University of Oregon high performance computer cluster, Talapas.
(for instance u → u + g). We denote the momentum and flavor splitting variables collectively by ζ. The function T l (ζ p , {p} m ) specifies our choice of the shower time t (Eq. (35)).
The last factor in Eq. (A1) contains color factors. The notation is from Ref. [4]. For instance, the operator t † l (f l →f l +f m+1 ) in t † l (f l →f l +f m+1 ) ⊗ t k (f k → f k +f m+1 ) acts on the ket color state and supplies the color matrices to split parton l into a new parton l and parton m + 1, with color representations according to the specified flavors. For a g → g + g splitting, we have (from Eq. (7.24) of Ref. [4]), where a † + (l) inserts the new gluon just to the right of parton l on whatever string in the color basis state contains parton l and a † − (l) inserts the new gluon just to the left of parton l. Here we define "left" to be the direction of the quark line in the color string. For the emission of a gluon from a quark line, we define (from Eq. (7.25) of Ref. [4]), where a † + (l) inserts the gluon at the quark end of the string to which the quark l belongs.
The functions λ lk ({p, f } m , ζ) are The function w dipole lk is the familiar eikonal splitting function, 3) of Ref. [7]. Here and throughout this paper, we assume massless partons.) The function A lk partitions the dipole splitting function, which is symmetric under k ↔ l, into a part considered to be a splitting of parton l and a part considered to be associated with a splitting of parton k. Our preferred choice is given in Eq. (7.12) of Ref. [6], (A6) The function w eikonal ll , from Eq. (2.10) of Ref. [5], is Here the total momentum of all the partons before the splitting is Q and after the splitting isQ. The tensor which is the polarization sum for gluon m + 1 in timelike axial gauge,Q · A(x) = 0. The function w ll is given in Refs. [4] and [5]. 10 It is rather complicated, but it is simple in the limit y → 0 with fixed z. For example, for a final state q → q + g splitting, we find (See Ref. [5], Eq. (2.24)) Here we recognize that (1 + z 2 )/(1 − z) is the DGLAP splitting function for this splitting. For a final state g → g + g splitting, we find (See Ref. [5], Eq. (2.52)) If we add this quantity and the same quantity with z ↔ (1 − z), corresponding to interchanging the two identical final state gluons with labels l and m + 1, the sum of the quantities in square brackets is the DGLAP splitting function for finding a gluon in a gluon. For an initial state q → q + g splitting, we find (See Ref. [5], Eq. (2.38)) Again, w ll is proportional to the DGLAP kernel (1 + z 2 )/(1 − z) for finding a quark in a quark. For an initial state g → g + g splitting, we find (See Ref. [5], Eq. (2.59)) (A12) Thus, the limiting form of w ll is proportional to the DGLAP kernel for finding a gluon in a gluon. Although w ll is proportional to the relevant DGLAP kernel in the limit y → 0 at fixed z, the full functions w ll are obtained directly from the relevant Feynman diagrams and are markedly different from the DGLAP kernels when y is not small, and particularly when y is comparable to z or 1 − z.
We have multiplied and divided by a factor N (k, l, ζ f ) that depends on whether k = l and on the flavors ζ f in the splitting: This factor, from Eq. (6.12) of Ref. [7], plays a role when we apply the LC+ approximation to the splitting.
In the case of gluon emission,f m+1 = g, the color factor G lc+ is non-zero only if parton k is color connected to parton l in the ket state or the bra state or both. That is, the result is zero unless the function χ(k, l, {c , c} m ) defined in Eq. (54) is nonzero.
The factor N (k, l, ζ f ) in Eq. (A13) is defined so that the total probability associated with the color factor G lc+ is just χ(k, l, {c , c} m ) times the probability associated with the starting color state: We have noted that ∆G(k, l, ζ f , {ĉ ,ĉ} m+1 , {c , c} m ) is non-zero only for k = l. We need also the operator V(t), which does not change the number of partons, their momenta, or their flavors and is related to H I (t) by 1 H I (t) = 1 V(t). Following Ref. [7] we define This creates a practical problem because the color operators are not diagonal in the basis that we use, so that it is not practical to calculate matrix elements of an exponential of V(t).
With the LC+ approximation, we define V lc+ (t), using 1 H lc+ (t) = 1 V lc+ (t). This gives us This corresponds to the way that we actually compute, but it is not suited for our present investigations. Instead, we replace × [t † l (f l →f l + g) ⊗ t k (f k →f k + g) + t † k (f l →f l + g) ⊗ t l (f l →f k + g)] .
The combination With a little manipulation, this becomes There is an overall factor 1/E 2 m+1 , so that there is a singularity in the soft limit. There is also an overall factor 1/( u m+1 − u l ) 2 , so that there would be a singularity in the collinear limit u m+1 → u l if this factor were not cancelled. The first term in h lk has a linear zero as u m+1 → u l , while the second term has a quadratic zero. This leaves an integrable collinear singularity, which disappears if one averages over the direction of ( u m+1 − u l ). Also, note that h lk can have either sign.
We could, if we liked, take the soft limit in w soft lk ({p,f } m+1 ) by settingp l → p l andp k → p k . For an initial state emission, we can also setQ → Q and we can replaceη → η in the parton distribution functions that multiply w soft lk ({p,f } m+1 ) in Eq. (B3). This doesn't help in the Deductor code, but it may be useful for analyses.