Parton shower evolution with subleading color

Parton shower Monte Carlo event generators in which the shower evolves from hard splittings to soft splittings generally use the leading color approximation, which is the leading term in an expansion in powers of $1/N_c^2$, where $N_c = 3$ is the number of colors. We introduce a more general approximation, the LC+ approximation, that includes some of the color suppressed contributions. There is a cost: each generated event comes with a weight. There is a benefit: at each splitting the leading soft$\times$collinear singularity and the leading collinear singularity are treated exactly with respect to color. In addition, an LC+ shower can start from a state of the color density matrix in which the bra state color and the ket state color do not match.


Introduction
Partons carry momentum, flavor, spin, and color. All of these quantum numbers are represented in parton shower Monte Carlo event generators like Pythia [1], Herwig [2], and Sherpa [3]. Spin and color do not fit easily into the event generator format because quantum interference between different spin and color states is important, even in the limit of very soft or collinear parton splittings. In this paper, we address the issue of color.
In previous work [4], we have shown how the evolution equations for a parton shower can be formulated in a way that fully includes spin and color. The resulting integrals can, in principle, be evaluated by Monte Carlo integration. However, simple Monte Carlo methods will not be practical when there are many splittings to be represented. In ref. [5], we found that these same evolution equations simplify if we average over spins and take the leading color limit, dropping all terms proportional to 1/N 2 c , where N c = 3 is the number of colors. Then one gets evolution equations that can be realized as a Markov process with positive probabilities. In ref. [6], we described a method for incorporating spin interference that we believe is practical. We are currently working on code to realize all of this. In order to keep the exposition in this paper as simple as possible, we average over spins and concentrate on color only.
The problem of treating color interference is, in our experience, more difficult. One method is to order gluon emissions according to emission angle, which takes into account quite a lot of the physics of color coherence at the cost of approximating the full dependence of the emissions on the emission angles [7]. With the approximations, this does give positive splitting probabilities. The other main method is to simply drop terms that correspond to interference in the color space on the grounds that these terms are suppressed by factors 1/N 2 c . This is the standard leading color (LC) approximation.
In this paper, we describe a less restrictive approximation, which we call the LC+ approximation. There is a cost to using the LC+ approximation in place of the LC approximation: one then gets shower events with weights, which can be negative. However, we argue that the deleterious effect of the weights on the numerical performance of the algorithm can be controlled -for instance by switching back to the LC approximation after some number of parton splittings. There are several advantages to the LC+ approximation. First, the color changes in collinear emissions and in emissions that are both soft and collinear are treated exactly. The only approximation is in wide angle soft emissions. Second, LC+ evolution can be started from any partonic color state, including states in which partons in the quantum amplitude have one color configuration while partons in the conjugate quantum amplitude have a different color configuration. The LC approximation cannot work with this kind of interference state. Third, the LC+ approximation is quite flexible, so that it could be applied within any parton shower program based on color dipoles, such as Pythia or Sherpa. For this reason, we believe that the method is of quite general interest.
We find that to state the method precisely, we need a fairly elaborate formalism based on linear operators on a certain vector space. We borrow this formalism from Refs. [4][5][6]. However, the essence of the LC+ approximation can be understood from a consideration of standard Feynman-like diagrams that represent color flow. For this reason, we provide a heuristic introduction to the approximation in sections 2, 3, and 4. We then turn to the more detailed, operator based, analysis in sections 5, 6, and 7. In the LC+ approach, there are different color states for the amplitude {c} m and the conjugate amplitude {c } m . In section 8, we describe how to get back to diagonal configurations, {c } m = {c} m at the end of the shower. In section 9 we describe how one can include in perturbation theory the effects that are neglected in the LC+ approximation. In section 10 we include the effects of the color phase induced by exchanging virtual soft gluons. Finally, in section 11, we conclude the main part of the exposition. We treat some more technical topics in three appendices.

Statement of the problem
In this section, we seek to describe how color evolves in a parton shower, illustrating why it is not simple to describe the color evolution exactly. First, let us note that, in a perturbative treatment, a parton carries momentum, flavor, spin and color. We can imagine that the momenta of partons at the end of the shower are quite precisely measured. Then the momenta of intermediate partons in a shower are also well determined. Thus we do not need to consider interference between states in which intermediate partons have different momenta p and p . Since the flavor of a mother parton in a splitting is determined by the flavors of its daughters, we also do not need to consider interference between states in which intermediate partons have different flavors f and f . However, this argument does not hold for spin and color. For instance, in a splitting 1 → 2 + 3 for gluons with colors a 1 , a 2 and a 3 , the matrix element is proportional to f a 1 ,a 2 ,a 3 and for fixed a 2 and a 3 , f a 1 ,a 2 ,a 3 can be nonzero for more than one value of a 1 . Thus we need to allow for quantum interference between different color states of the same parton. We also need to allow for quantum interference between different spin states.

Evolution without color or spin
Let us consider what happens in a parton shower that evolves from hard splittings to soft splittings. 1 To get started, we ignore both spin and color. We define a "shower time" variable t such that an initial hard parton scattering happens at t = 0 and then at each interval dt a parton has some probability of splitting to become two partons. Harder splittings happen at smaller t values, successively softer splittings happen at larger t values. For instance, exp(−t) can be proportional to the virtuality in a splitting or to the transverse momentum of the daughters relative to the mother parton direction. As implemented in a computer program, the partonic system always has a definite state. Ignoring spin and color, the state of the system is described by the momenta and flavors, Here p a is the momentum of the incoming parton from hadron A, equal to a fraction η a of the hadron momentum, p b is the momentum of the incoming parton from hadron B, and the p i are the momenta of m outgoing partons. The flavors of the partons are denoted by discrete flavor variables f . In each shower time interval dt, there is a certain probability that the system will switch to a new state. What actually happens in a given simulated event is determined by generated pseudo-random numbers. This means that in an ensemble of simulated events, there is an probability ρ({p, f } m , t) that the partonic system is in a certain state at shower time t. This probability distribution then evolves with t. Here H I (t) is a linear operator on the space of probability distributions. The operator H I (t) corresponds to the splitting probabilities chosen for shower evolution. (There are many possible choices.) Then V(t) is another linear operator that is constructed from H I (t) so as to conserve probability in the shower evolution. The Sudakov factor that represents the probability that there was no splitting between times t 1 and t 2 is exp(− t 2 t 1 dt V(t)). We will return to this in later sections. For the moment, all that we need to know is that the evolution of the probability function ρ is determined by parton splitting probabilities.

Ignoring spin
Now we need to consider spin and color. We have addressed the spin problem in ref. [6]. Combining the spin treatment of ref. [6] with the discussion of this paper is not difficult, but adds a layer of complexity. In order to keep this paper as simple as possible, we ignore spin by supposing that we sum over the spins of the daughter partons in a splitting and average over the spins of the mother parton. Thus in this paper we address the color problem but not the spin problem.

The color density operator
The natural language for describing an evolving probability distribution is statistical mechanics. Indeed, eq. (2.2) is a standard sort of equation for evolution in statistical mechanics. In this paper, we want to include quantum color, so we need quantum statistical mechanics. Thus we need a density operator that depends on t and is an operator on the  Our analysis in this paper makes use of the standard color basis defined by two sorts of vectors [4,8]. First, there are open string vectors Ψ(S) {a} = N (S) −1/2 [t a 2 t a 3 · · · t a n−1 ] a 1 ,an . (2.4) Here a 1 is a quark color index, a n is an antiquark color index, and the other a i are gluon color indices. The t a are standard SU(3) color matrices in the fundamental representation. Also, N (S) = N c C n−2 F is a normalization factor such that Ψ Ψ = 1. Second, there are closed string vectors Here all of the a i are gluon color indices. Also, N (S) = C n F is a normalization factor such that Ψ Ψ = 1 − [−1/(N 2 c − 1)] n−1 . This normalization factor is close to 1, with a small correction. The most general color basis vector, which we denote by {c} m , is a product of these two kinds . That is, the basis vectors are orthogonal in the N c → ∞ limit. Now let's look at an example. In figure 2, we depict {c} m {c } m for a color structure that arises in the hard scattering processq(a)+q(b) → g(1)+g (2). Both {c} m and {c } m show the color state that corresponds to t-channelq exchange. Thus we have a diagonal contribution to the color density operator, with {c} m = {c } m .   One can also have u-channelq exchange, which amounts to exchanging gluons 1 and 2. In figure 3, we still have the t-channel diagram for {c} m , but now we illustrate the contribution from the u-channel diagram for {c } m . This gives an off-diagonal contribution to the color density operator, with {c} m = {c } m . Each of the two contributions to the color density operator shown in the two figures can serve as the starting point for a parton shower. Since their color structures are different, they should generate different showering.

Color structure of the parton splitting operator
Having seen the meaning of the color density operator, we can now consider what happens in shower evolution when a gluon splits. In figure 4, we show the color structure when mother gluon 1 splits into daughter gluon 1 and daughter gluon 3. Using the identity if a 1 a 3 c t c = t a 1 t a 3 − t a 3 t a 1 , we find that there is a term in which the new gluon 3 is inserted below gluon 1 on the quark line minus a term in which the new gluon is inserted above gluon 1 on the quark line. 2 Applying this identity, we see in figure 5 that letting mother gluon 1 split in both the ket state and the bra state in figure 2 gives contributions with four new {c} m+1 {c } m+1 terms. Even when we use only the simplest kind of splitting and even though we start with a color diagonal density operator, we generate off diagonal terms in the color density operator.
There can also be quantum interference between emissions of a soft gluon from two different partons. We illustrate this in figure 6. Mother parton 1 emits gluon 3 in the ket state, while the antiquark with label b emits gluon 3 in the bra state. In a physical gauge, the matrix elements for this are large as long as gluon 3 is soft. One can consider that a color dipole consisting of partons 1 and b emits the gluon. In a dipole shower, one wants to account for emission from each dipole (insofar as complications from color allow). In a dipole antenna shower like Vincia [9], one keeps the dipoles together as a unit. In this paper, we follow the method of a partitioned dipole shower like Pythia, which partitions the radiation pattern into two terms. One term is singular when the soft gluon is collinear with parton 1, the emitting parton, but is not singular when the soft gluon is collinear with parton b, the helper parton. The other term is singular when the soft gluon is collinear with parton b but is not singular when the soft gluon is collinear with parton 1. For our present discussion, let us consider parton 1 to be the emitting parton. On the right hand side of the figure, we have applied the color identity of figure 4 to expand this contribution in our standard color basis. We see that we get two off diagonal contributions. Parton 3 can also be emitted from parton 1 with parton 2 as the helper parton. This case is illustrated in figure 7. When we expand in five-parton color basis states, there are four contributions. One of them is color diagonal.
Let us summarize. Using a few simple examples, we have explored the color content of the splitting operator H I (t) in eq. (2.2). There are two sorts of graph, direct graphs like figure 5 and interference graphs like figures 6 and 7. The color structure of both kinds of graphs is non-trivial. Even if we start from a color diagonal contribution to the color density operator, each splitting creates off-diagonal contributions. With each new splitting, the color state becomes more complicated.

Color structure of the virtual splitting operator
What about the operator V(t) in eq. (2.2)? This operator maps a state of the color density operator with m final state partons into another state with m final state partons. Physically, it represents virtual Feynman graphs with hardness scale characteristic of the current shower time t. We determine what V(t) should be by demanding that shower evolution conserve probability. This means that where multiplying by 1 is a convenient way of writing the instruction to take the trace of the color statistical operator and integrate over all of the parton variables. Note that taking the trace of the color statistical operator means using (2.8) This notation denotes that (h + iφ) is an operator on the ket part of the color density operator and (h † − iφ) is an operator on the bra part. The operator φ is hermitian and represents a color phase. We will consider φ in section 10. Until then, we simply set φ = 0. The operator h is hermitian in the full color theory, but the LC+ approximation for it is not. For that reason, we distinguish the roles of h and h † in our formulas. With φ = 0, we have We can state this a little more precisely. If ρ({p, f } m , t) is defined by eq. (2.3) and ρ = V ρ , then (2.10) Each term in H I determines a contribution to h. To see how this works, consider the splitting in figure 5. We define h so that the corresponding contribution to ρ has the color factor shown in figure 8. To verify probability conservation, eq. (2.6), just take the color trace of figure 8 and compare it to the color trace of the left hand side of figure 5. In each case we get the number represented by the color diagram in figure 9.
The color structure of V is more complicated in the case of interference diagrams. Consider the splitting in figure 7, in which parton 1 in the ket state splits and we take the   interference with the emission of the same gluon from helper parton 2 in the bra state. We can consider this diagram to correspond to a contribution to V in which h acts on the bra state. We define this contribution to h so that ρ has the color factor shown in figure 10. To verify probability conservation, just take the color trace of figure 10 and compare it to the color trace of the left hand side of figure 7. Then the contribution to V in which h acts on the ket state corresponds to the complex conjugate of the splitting in figure 7.
Note that the contribution to V(t) corresponding to a direct splitting term in H I is trivial. Evidently from figure 8, this contribution to V(t) simply multiplies the color state by C A . In contrast, the contribution to V(t) corresponding to an interference term in H I is not trivial. When one expands the bra state in figure 10 in the color basis states, we get several contributions. 3 Thus, in general, V(t) is not diagonal in color. This makes an exact accounting for color difficult because V(t) enters the generation of parton showers in the form of the Sudakov factor, the time ordered exponential T exp −

The leading color approximation
We have now seen something of the structure in color space of a leading order parton shower (of the partitioned dipole variety). Of course, this is not the structure of any existing computer program that models a parton shower with quantum color. There is no such program. The evolution equations make good sense, but implementing more than a couple of splitting steps as a computation is not feasible with current methods. There is, however, a simple approximation that one can apply to get a practical implementation. This is the leading color (LC) approximation.
The general idea of the leading color approximation is to neglect contributions to the probability to get a given final partonic state that are suppressed by powers of 1/N 2 c . To proceed, one replaces gluons in the 8 representation of SU (3) by gluons in the 3 ×3 representation, using  figure 5, one keeps the first two terms and drops the second two terms. On the right hand side of figure 7, one keeps the first term and drops the other three. On the right hand side of figure 6, one drops both terms. That is, interference between emission of a gluon from parton l and from a helper parton k that is not directly color connected to l is neglected.
With the leading color approximation, the no-splitting operator V is simple. For instance, on the right hand side of figure 5 we keep the first two terms. In each of these terms, when we take the color trace we see that adding the emitted gluon 3 creates one more color loop and thus one more factor of N c . Thus we can take the color operator in V to simply multiply the color state we started with by C A = N c (including two graphs, each with a factor 1/2 from eq. (3.1)).
For a splitting q → q + g, the color factor in V would be N c /2 when we calculate this way, but one normally uses C F instead by simply multiplying the splitting probability by Strictly speaking, one would drop all g → q +q splittings in the leading color approximation, but one normally keeps the part of this splitting suppressed by 1/N c , omitting the parts suppressed by 1/N 3 c . The leading color approximation gives a simple shower algorithm. It is also intuitively appealing. The ket states and the bra states in the color density operator always have the same {c} m . In {c} m , we can think of each quark or antiquark as being connected by a color string to a neighboring gluon in a color basis state in figure 1. Each gluon is connected to two other partons by a color string. Then emitting a new gluon means connecting it to two of the previous strings. Interference between emitting a new gluon from parton l and helper parton k can only occur if the two partons were color connected; then the new gluon is connected to the string that previously joined partons l and k.
One can note two problems with the leading color approximation. First, it neglects terms suppressed by powers of 1/N 2 c . Second, it cannot start with color density operator contributions {c} m {c } m with {c} m = {c } m .

Introduction to the LC+ approximation
We propose in this section an improved "LC+" approximation that goes beyond the leading color approximation by including some of the contributions to cross sections that are suppressed by powers of 1/N 2 c . To go beyond the leading color approximation, one has to give up something. 4 We choose to give up having an algorithm that can be implemented without having weights for events. In particular, the terms in the splitting probabilities in the LC+ approximation have both plus and minus signs. One cannot generate events with negative probabilities, so in order to include these terms one will have to generate certain splittings with positive probabilities and negative weights. The weights are then carried with the event. Now, having weighted events is not intrinsically a problem for convergence of Monte Carlo integrations to a physical answer. However, there would be a problem that would not be easy to avoid if a physical cross section to be calculated had the form σ 0 [1.00 − 9.01/N 2 c + · · · ]. Then calculating the result by Monte Carlo integration would not give an accurate answer for N c = 3. We believe that this does not happen. If it does happen for some observable cross section, then the LC+ approximation may not work well. Of course, in this case, the standard LC approximation is giving us the wrong answer and the LC+ result would at least give an indication of problems.
We will state what the LC+ approximation is quite precisely using an operator notation in section 6. However the main idea can be grasped quite easily using the examples that we have already explored, so we do that in this section.

The parton splitting operator
First, the direct graph for the splitting of parton l by emitting a gluon labeled m + 1 contains (in a physical gauge) the collinear singularity that occurs when the parton momenta after the splitting obeyp m+1 ∝p l and also the double singularity that occurs when p m+1 → 0 withp m+1 ∝p l . We want the coefficients of the logarithms that arise from these singularities to be exactly right with respect to color. Thus we keep the color structure of these splittings exactly. This means that in the right hand side of figure 5 we keep all four terms. We can start with {c } m = {c} m as in figure 11. Again, there are four terms and we keep them all.  Now consider a splitting in which gluon emission from parton l in the ket state interferes with emission from helper parton k in the bra state. An example was shown in figure 6. Here the helper parton (k = b) is not color connected to the emitting parton (l = 1) in the bra state. In this case, the LC+ approximation is to drop the contribution entirely. A second example was shown in figure 7. Here the helper parton (k = 2) is color connected to the emitting parton (l = 1) in the bra state. In this case, the LC+ approximation is to keep the contributions in which the emitted gluon attaches between l and k in the bra state. That is, we keep the first and third terms on the right hand side of figure 7. A third example is shown in figure 12. In this case, {c } m = {c} m . The helper parton (k = 2) is color connected to the emitting parton (l = 1) in the bra state, so the LC+ approximation is to keep the contributions in which the emitted gluon attaches between l and k in the bra state. The two contributions that are retained in the LC+ approximation are shown on the right hand side of figure 12.

The virtual splitting operator
To find the structure of the virtual splitting operator V(t) in the LC+ approximation we use the definition (2.9). For a direct splitting diagram, we consider the real emission diagram to correspond to the sum of virtual diagrams in which h acts on the ket state and in which h † acts on the bra state. For an interference diagram with the helper parton in the bra state, we consider the real emission diagram to correspond to a virtual diagram in which h acts on the bra state. An interference diagram with the helper parton in the ket state corresponds to a virtual diagram in which h acts on the ket state. Then we impose probability conservation, eq. (2.6), to relate V(t) to the definitions of the LC+ approximation for H I (t) outlined above.
For a direct splitting diagram as in figure 11, this gives the contribution to V(t) illustrated in figure 13. To see this, we simply note that the color trace of the left hand side of figure 11 is the same as the color trace of figure 13. For the interference diagram in figure 12, we first note that the two terms that are retained in the LC+ approximation sum to the contribution shown in figure 14. This gives the contribution to V(t) illustrated in figure 15 since the color trace of the right hand side of figure 14 is the same as the color trace of figure 15.
We now note something remarkable. The color basis states are eigenstates of the operator h that defines V(t). For a direct splitting of a gluon, as in figure 13, the eigenvalue is C A /2. Similarly, for a direct splitting of a quark or antiquark, the eigenvalue is C F /2. For an interference diagram in which a gluon splits to two gluons with interference from gluon emission some other parton, as in figure 15, a simple calculation shows that the eigenvalue is C A /4. Similarly, for an interference diagram in which a quark or antiquark splits, emitting a gluon with interference from gluon emission some other parton, the LC+ approximation is shown in figure 16. Then the corresponding contribution to V(t) has eigenvalue C F /2.

Operator based analysis with full color
In the preceding sections, we have sketched the main idea of the LC+ approximation. Now we need to make the idea precise. To do that, we first set up a formalism to describe the evolution of generic parton shower that includes quantum color exactly. To keep the presentation simple, we average over spins. More precisely, we sum over the spins of partons after each splitting and average over the spins of the partons that enter each splitting. Then spin is not visible in the evolution equations at all. The formalism is taken from Refs. [4][5][6]. However, here we keep the shower quite generic by not specifying exactly the splitting functions, the definition of the shower time t that is used to order successive splittings, or the momentum mapping for each splitting that allows all partons to be on-shell while at the same time exactly conserving momentum.

Parton labels
At each stage of the shower, there are two initial state partons with labels "a" and "b" together with m final state partons with labels 1, . . . , m. The parton momenta and flavors are then specified by a list {p, f } m as in eq. (2.1). 5 At each step of the parton shower, any one of the partons can split. This includes an initial state parton, which splits in "backwards evolution" to another initial state parton plus a final state parton. In either case, let l be the label of the parton that splits. At a splitting, the parton l remains and one more parton, with label m + 1 is created. After the splitting, the momenta and flavors are {p,f } m+1 . Whenever a final state gluon is created, we assign the label m + 1 to the gluon. In the case of a final state g → g + g splitting, one daughter gluon has label l and the other has label m + 1. We use the interchange symmetry of the process to rearrange the splitting function so that there is a singularity when gluon m + 1 becomes soft but no singularity when gluon l becomes soft.
In the case of gluon emission, there are interference graphs. A gluon m + 1 emitted from parton l in the ket state can be emitted by partner parton k in the bra state. Similarly a gluon m + 1 emitted from parton l in the bra state can be emitted by partner parton k in the ket state. The interference diagrams are important when gluon m + 1 is soft.

The evolution equation
The development of a parton shower can be described by giving an evolution equation for the operator ρ({p, f } m , t) defined by eq. (2.3) as the density operator on the color space of the evolving partonic system. (It would be an operator on the spin space also except that we average over spins in this paper.) The density operator is a function of the parton momenta and flavors. This function can be regarded as a vector ρ(t) in the space of functions of {p, f } m with values in the space of operators on the quantum color space. We take the vector ρ(t) to obey a linear evolution equation of the form (2.2).
In eq. Making use of these basis vectors, we can define H I (t) by specifying the matrix elements of H I (t) between a partonic state after a splitting and the state before the splitting, {p,f ,ĉ ,ĉ} m+1 H I (t) {p, f, c , c} m . We begin this task in the next few subsections by discussing some of the ingredients in this matrix element.

Splitting functions
To describe emission of a soft gluon m + 1 from parton l with interference from emission from helper parton k, we use the dipole splitting function In this expression, partons l and k can have nonzero masses. The expression for w lk may be more familiar in the massless limit,p 2 l → 0 andp 2 k → 0, where it becomes Eq. (5.2) includes all four diagrams for emission from either parton l or parton k in both the bra and ket states, calculated in the limitp m+1 = λP m+1 with λ → 0. If we calculate the four contributions separately using the eikonal approximation in a physical gauge, then The first two terms represent the direct graphs while the last term represents the two interference graphs. 6 The splitting function w eikonal ll in eq. (5.4) is singular in the limit in which parton l is massless and parton m + 1 is collinear with parton l:p m+1 → [(1 − z)/z]p l . In this limit, There is a factor of the virtuality of the splitting, 2p m+1 ·p l in the denominator and there is a function of the momentum fraction z that is singular in the soft gluon limit z → 1. The behavior of this function away from the collinear limit depends on the conventions used to define the shower. For instance, it depends on the definition of z.
The splitting of parton l into a parton of the same flavor plus a gluon both in the ket state and in the bra state is described by a function w ll . Suppose that parton l is a quark. Then in the limit in which parton l is massless and parton m + 1 is collinear with parton l we have Away from this limit, the exact result depends on the conventions used to define the shower. 7 Note that w eikonal ll matches w ll in the limit of a collinear splitting that is also soft, For other sorts of splittings, there are splitting functions w ll that are proportional to the DGLAP splitting kernels in the limit of massless, collinear splittings.

Shower time
Parton showers are based on evolution of the system as a shower time variable t increases. The idea is to start at the hard interaction and move to a first splitting that is less hard, then move to softer splittings as t increases. For initial state splittings, this means moving backward in physical time.
As a measure of softness, one can use the virtuality of the splitting, the virtuality divided by the energy of the mother parton, or the square of the part of the momentum of one of the daughters that is transverse to the direction of the mother. The shower 6 The function w eikonal lk , which equals w eikonal kl , is denoted by W lk /(2A lk ) and given in eq. (4.4) of ref. [5]. The function w eikonal ll is denoted by W eikonal ll and given in eq. (2.10) of ref. [5]. These functions include only the leading λ −2 contribution in the limit λ → 0. One could, in principle, add contributions proportional to λ −1 and higher powers of λ. 7 For a g → g + g splitting, the DGLAP splitting kernel appears in w ll , as in eq. (5.6), once one symmetrizes in z ↔ 1 − z. Before symmetrization, the result depends on the conventions used to define the shower.
program Herwig orders splittings according to the angle between the daughters times the energy of the mother, with wide angle splittings first. This ordering has its advantages, but the generic shower scheme outlined here would need some modification to fit with angular ordering.

Momentum mapping
At each splitting, one starts with m final state partons with momenta {p} m and ends with m + 1 final state partons with momenta {p} m+1 . One might like to definep j = p j for j = l and j = m + 1 and to takep l +p m+1 = p l (orp l −p m+1 = p l for an initial state splitting). However, this would not allow all three ofp l ,p m+1 , and p l to be on-shell. Accordingly, one needs to take some momentum from the partons other than l in order to supply the needed momentum in the splitting and keep exact momentum conservation. Thus we need a momentum mapping in which the {p} m plus three splitting variables determine the {p} m+1 . The three splitting variables can be, for instance, the virtuality (p l ±p m+1 ) 2 − m 2 (f l ), a momentum fraction z and an azimuthal angle of the daughters about the direction of the mother. Equivalently, the {p} m+1 determine three splitting variables and the {p} m . The needed momentum mapping can be specified by giving a function {p,f } m+1 P l {p, f } m that consists of an integration over the three chosen splitting variables and a product of delta functions that determine the {p} m+1 from the splitting variables and the {p} m . There are many ways to do this, but for our purposes, all we need to know is that defining a shower evolution entails the specification of For our generic shower, we take the mapping operator P l to depend on the label l of the parton that splits. One can also let the mapping operator depend on the index k of the helping parton involved in interference diagrams. This is a common choice because it allows all of the momentum transfer to come from changing the momentum of parton k. That is, one can takep l +p m+1 − p l = p k −p k . However, this scheme is a bit awkward for cases in which there is no helper parton k, such as g → q +q. Accordingly, we keep the generic shower equations simple by letting P l depend only on the index l.

Splitting operator
Now we are ready to specify the matrix elements of H I (t): In the first line on the right hand side of eq. (5.7), we have a sum over indices l and k of partons. The parton with label l is the one that splits. There is another label k so that we can include graphs that represent quantum interference between emission of a gluon from parton l and from another parton k. We call parton k the helper parton. For the quantum interference terms, we have k = l. There are also graphs that do not represent quantum interference. For these, k = l.
The first line on the right hand side of eq. (5.7) also includes a delta function that specifies the definition of the shower time t, then the function that defines the momentum mapping.
In the second line, we have a ratio of parton distribution functions. For a final state splitting, this ratio is 1. For an initial state splitting, this ratio replaces the parton distribution functions at the previous momentum fraction η a or η b by the parton distribution function at the new momentum fraction after the splitting.
In the following lines of eq. (5.7) there are three terms, corresponding to three types of splittings.
First, there are splittings in which parton m + 1 is not a gluon. For example, a final state g → q +q splitting is in this class. For these cases, there is a splitting function w ll ({p,f } m+1 ) that is singular for a collinear massless splitting, but does not have a soft parton singularity.
In the second term in eq. (5.7), there are splittings in which parton l emits a gluon. The splitting function here, w ll − w eikonal ll , is singular for a collinear splitting, but the singularity when the emitted gluon is soft has been removed by the subtraction.
The most interesting term in eq. (5.7) is the third, in which a soft gluon m+1 is emitted from a dipole consisting of partons l and k. The splitting function w dipole lk is calculated using the eikonal approximation according to eq. (5.2). It is singular when the gluon is soft and also contains singularities when the direction of the gluon momentum is collinear to either the direction ofp l (ifp 2 l = 0) or the direction ofp k (ifp 2 k = 0). It is symmetric under interchange ofp l andp k .
We have manipulated the third term in order to separate the roles of partons k and l. Let A lk be a function of the momenta {p} m+1 and A kl be the same function with the roles ofp k andp l reversed. Furthermore, let A lk + A kl = 1. Then, by interchanging the names of the dummy indices l and k, we have This gives the factor A lk that multiplies w dipole lk in eq. (5.7). If we were to take A lk = 1/2, this manipulation would do nothing, but instead we choose A lk so that A lk = 1 when p m+1 is parallel top l and A lk = 0 whenp m+1 is parallel top k . Then A lk w dipole lk is singular when the gluon is soft or collinear with parton l but not when it is collinear with parton k. Our preferred choice for A lk is given in eq. (7.12) of ref. [6]. One can also use the Catani-Seymour dipole partitioning function [10] for this purpose.
Finally in eq. (5.7) there is a color factor, which is the main focus of this paper. There are two terms in the color factor, which are related by interchanging the indices l and k. If k = l, the two terms are identical. In the first term, there is an operator t † l (f l →f l +f m+1 ) that acts on color ket states {c} m . This operator attaches the proper color matrix for the splitting to the color index for parton l. Similarly, we apply the operator t k (f k →f k +f m+1 ) to the bra color state {c } m . The resulting color bras and kets can be expanded in color basis states: Then the matrix element that we want is the corresponding expansion coefficient:

The virtual splitting operator
The virtual splitting operator V(t) leaves the number of partons unchanged and does not change their momenta or flavors. It does, however, multiply by a matrix in color space. This matrix has two terms. First, there can be a correction to the ket color state with no change to the bra colors. Then, there can be a correction to the bra color state with no change to the ket colors. Thus the color structure of V(t) is the color structure of one loop virtual corrections. We write The color matrices H L and H R are the matrices that represent operators h + iφ and [h + iφ] † = h † − iφ that act on the ket color vectors and the bra color vectors, respectively: Here φ = φ † . In the simplest formulation, φ = 0. We consider φ = 0 in section 10. In the full color treatment, we will have h † = h. However in the LC+ approximation we will It is useful to solve the shower evolution equation (2.2) in the form Here N (t 2 , t 1 ) is the no-splitting operator, that provides the Sudakov factor that is usually interpreted as the probability to not have a splitting between time t 1 and time t 2 . In N (t 2 , t 1 ), the virtual splitting operators V(τ ) are time ordered. Iterating eq. (5.14) gives ρ(t) as a power series in H I , with H I evaluated at times t i and factors of N (t i+1 , t i ) in between. We interpret eq. (5.14) as saying that possibly the system gets from t 0 to t without splitting. If not, it goes from t 0 to τ without splitting, then splits according to H I (τ ), then evolves further according to the full evolution operator U(t, τ ). Given the structure (5.11) and (5.12) of V(τ ), the operator N (t i+1 , t i ) has the form of a matrix in labels of the color basis elements, Here the matrix elements are obtained from the operators h that define V(t), Thus the Sudakov factor breaks into two factors, one for the ket state and one for the bra state. That is, there is Sudakov factor for each color ordered quantum amplitude. That is remarkably simple. However, the structure of the Sudakov factor for each quantum amplitude is remarkably complicated since it is an operator that mixes color basis states.

Probability conservation
We relate V(t) to H I (t) using the requirement that showering not change the probability of the hard scattering event that initiates the shower. Probability conservation requires eq. where, to be precise, We are particularly interested in the color structure of this. We note by examining eq. (5.7) that on the right hand side of this equation the following generic color factors will occur: If we take the color trace of eq. (5.9) and compare it to the trace of eq. (5.10) we have If we use eq. (5.7) in eq. (5.23) and use the identity (5.24), we find (5.25)

Structure of the virtual splitting operator
Given the structure of V(t) in eq. (5.11) and the definition (5.19) of 1 , the right hand side of eq. (5.18) is Here the operator on the right hand side is (h + iφ) + (h † − iφ) = h + h † . Thus this relation determines h + h † but not φ. We have assumed that φ = 0. In section 10, we will explore the color phase φ. We demand probability conservation, eq. (5.18). Comparing eq. (5.27) to eq. (5.25), we see that probability conservation holds if In the color factor, we have a factor t l t † k . Then h † is given by the same expression with a factor t k t † l . In fact, h = h † because t k t † l = t l t † k . For k = l, this is obvious. For k = l, we are dealing with a gluon with color index a exchanged between the lines k and l. We note that t k t † l inserts SU (3) generator matrices T a in the appropriate representations on lines l and k, then sums over a. The generator matrices are self adjoint and they commute with each other because they act on different parton lines, so that t k t † l = t l t † k . This result might perhaps be regarded as elegant, but it poses difficulties for practical implementation in a parton shower Monte Carlo program. The difficulty comes from the fact that the operator t l (f l → f l +g) t † k (f k → f k +g) is represented by a non-diagonal matrix in the standard color basis. In the full h({p, f } m ) we have a sum of such operators with momentum dependent coefficients. For the Sudakov factor (5.17), we need the exponential of an integral over t of this matrix. Thus the Sudakov factors are complicated and it is not easy to see what to do with them.

Color suppression index
To help with our analysis, we define a quantity P that we can call the color suppression power and a related but more useful quantity I that we can call the color suppression index.
To define the color suppression power, we first note that shower evolution can generate a factor 1/N pc c that comes from directly from applying the splitting operator to the color states {c} m and expanding in color basis states {ĉ} m+1 . Thus p c is the power of 1/N c that is associated with the quantum amplitude.
For all of this paper except section 10, the sole source of non-zero p c is g → q +q splittings. In a g → q +q splitting, there are two ways to connect the new q andq to the previous color state. The original gluon color was defined by inserting t a ij along a color 3 line. Now with the g → q +q splitting we have a color matrix t a ij t a i j in the color amplitude. We can use the Fierz identity, where c P (m) = 0. The color suppression power is of obvious interest, but is less useful than we would like because cancellation among terms in the expansion of {c } m {c} m can lead to P being larger than it is for individual terms in the expansion.
We can improve on the definition by defining a color suppression index I according to That is, if we use I(m) to estimate the amount of color suppression, we can never overestimate. Second, at each stage of the shower, the color suppression index either stays the same or else it increases: Thus we can think of I as measuring color disorder, like entropy: it can never decrease as the shower evolves. Both of eqs. (5.32) and (5.33) can be proved in a fairly straightforward way. We omit the proofs.

How the color suppression index grows
In this section, we investigate how the color suppression index I(m) grows. Imagine calculating {c } m {c} m U (Nc) . Look at any one of the gluons in the state, say gluon l. Its color factor was a t a i i t a j j . With the U (N c ) approximation, this becomes (1/2) δ i j δ j i . Counting the rest of the graph, the indices can be connected in two possible ways In the first case, let us call gluon l a healthy gluon. In the second case, let us call gluon l a frail gluon. 8 With a healthy gluon, we get a factor N 2 c for the sum in eq. (5.34). With a frail gluon, we get a factor N 1 c for the sum in eq. (5.35). Let us look at this in more detail. Consider eq. (5.34) for the healthy gluon. If the healthy gluon were not there, we would have a factor δ i i δ ii = N c . Now insert the healthy gluon and divide by C F ∝ N c to normalize the states with one more gluon according to eqs. (2.4) and (2.5). The new overlap then has a factor N 2 c /N c = N c . That is, the color suppression index has not changed. Now consider eq. (5.35) for the frail gluon. If the frail gluon were not there, we would have a factor δ ii δ jj = N 2 c . Now insert the frail gluon and divide by C F ∝ N c to normalize the states with one more gluon. The new overlap then has a factor N c /N c = 1. That is, the color suppression index has increased by 2.
We can now look at what happens when a frail gluon splits and what happens when a healthy gluon splits.
Suppose first that gluon l is a frail gluon and splits, creating another gluon m + 1. Within the LC+ approximation, there are two configurations for the two gluons. The first is Suppose next that gluon l is a healthy gluon and splits. Within the LC+ approximation, there are again two configurations for the two gluons. The first is the parallel configuration, a,b [t a t b ] i i [t b t a ] j j . With this configuration, an easy calculation shows that both gluons in the new state are healthy gluons and the color suppression index is unchanged: I(m + 1) = I(m). The two gluons can also be in the crossed configuration, With this configuration, we find that both gluons in the new state are now frail gluons and the color suppression index of the state increases: I(m+1) = I(m)+2.
When gluon l splits, creating another gluon m + 1, there will normally be other gluons already present and the splitting of gluon l can change the health status of some of these other gluons. Consider another gluon with label s. When gluon l splits to gluons l and m + 1 in the parallel configuration, so that I(m + 1) = I(m), a simple argument shows that the health status of gluon s remains the same. Now suppose that gluon l is a frail gluon and splits to gluons l and m + 1 in the crossed configuration, so that I(m + 1) = I(m).
Then if s was healthy it remains healthy and if s was frail it may remain frail or it may become healthy, depending on how it was connected. Finally, suppose that gluon l is a healthy gluon and splits to gluons l and m + 1 in the crossed configuration, so that I(m + 1) = I(m) + 2. Then if s was frail it remains frail and if s was healthy it may remain healthy or it may become frail, depending on how it was connected. In summary, then, when a new gluon m + 1 is added so that it is color connected to gluon l in both the ket state and in the bra state, then splittings that increase the color suppression index may change healthy gluons to frail gluons while splittings that leave the color suppression index unchanged may change frail gluons to healthy gluons.
We will use these observations in appendices A and B.

The LC+ approximation
With the generic structure of a (spin averaged) leading order shower set up, it is pretty simple to define the LC+ approximation. In the definition (5.7) of H I (t), in the terms that represent quantum interference between emitting a gluon from parton l and emitting the same gluon from helper parton k we replace Here C(l, m+1) acting on a state {ĉ} m+1 gives 1 if partons l and m+1 are color connected in {c} m and 0 otherwise. Thus in the LC+ approximation we keep only the color states in which partons l and m + 1 after the splitting are color connected. We generalize this to include the possibility that k = l by replacing with a generalized definition of C(l, m + 1). For k = l, whenever parton m + 1 is a gluon, partons m + 1 and l after the splitting are always color connected, so we want C(l, m + 1) to return the color state unchanged. Also, for k = l wheneverf l = g, partons m + 1 and l after the splitting are also color connected, so we want C(l, m + 1) to return the color state unchanged. 9 In the case k = l withf l = q andf m+1 =q, or vice versa, the two daughter partons may not be color connected. We simply define C(l, m + 1) to give one in this case. Thus we define The operator C(i, j) is a projection operator but it is not an orthogonal projection operator because the basis states are not orthogonal. Thus C(i, j) = C(i, j) † .
Thus the LC+ approximation for H I is only slightly modified from the full H I in eq. (5.7): Here M c is the color matrix defined by The effect of the projection operator C † (l, m + 1) was illustrated in figure 7. With full color, the interference graph shown gives four contributions {ĉ} m+1 {ĉ } m+1 , but the projection operator C † (l, m + 1) removes the second and fourth contributions, in which gluon m + 1 = 3 is not color connected to parton l = 1 in the bra state. Only the first and third contributions remain. Similarly in figure 6 there are two contributions with full color but both are eliminated in the LC+ approximation. There are two terms in eq. (6.4) with k = l, one withf m+1 = g and one withf m+1 = g. The corresponding splitting functions have collinear singularities but not soft singularities. In both these terms, the projection operator C(l, m + 1) acts as the unit operator, so no approximation is made in these terms.
In the k = l term,f m+1 = g and the effective splitting function A lk w dipole lk is singular when the gluon is soft, collinear with parton l, or both. In the limit that the gluon is collinear with parton l or both collinear and soft, the LC+ approximation becomes exact, even though these terms contain the operator C(l, m + 1). To see this, note that A lk w dipole lk is independent ofp k in the limitp m+1 → λp l . Thus all of the terms with different indices k have the same coefficient A lk w dipole lk . Because of that, we can use the color identity l t l = 0 to write However, after a gluon emission from line l, parton m + 1 is always color connected to parton m + 1 in the new color state. Thus That is, the operator C has no effect in the collinear or soft×collinear limit. We conclude that the LC+ approximation becomes exact in the collinear or soft×collinear limit. It is only an approximation when an emitted gluon is soft but not collinear to the emitting parton.
For the virtual splitting operator, we use the the definition specified in eqs.
Now the operator h LC+ contains operators that act on the color space. However, the color basis vectors are eigenfunctions of these operators. In the case k = l, for which C(l, m + 1) acts as the unit operator, the operator t l (f l →f l +f m+1 ) t † l (f l →f l +f m+1 ) has eigenvalue T R , C F , or C A depending on the flavors of f l andf m+1 andf l . (The C A case is illustrated in figure 13.) For the case k = l, the operator t l (f l → f l + g) C(l, m + 1) t † k (f k → f k + g) is a little more complicated. Because of the projection operator C(l, m + 1), this operator gives zero unless the helper parton k is color connected to the emitting parton l in the state {c} m on which h LC+ acts. When l and k are color connected, one gets an eigenvalue C A /2 or C F depending on the flavor of parton l. (The C A /2 case is illustrated for h † in figure 15.) Thus where the eigenvalue λ LC+ is When l and k are color connected, the eigenvalue depends on whether k = l and on the flavors in the splitting: Since the color basis vectors are eigenvectors of h LC+ , the Sudakov factors for the color ordered amplitudes, eq. (5.17), are simply numerical factors: This is just the same as in the leading color approximation. In fact, with a suitable adjustment of what numerical factors one uses in the leading color approximation, the LC+ Sudakov factors are exactly the square root of the Sudakov factor for the leading color approximation.

Weights
Let us now look at a splitting step in the evolution equation (5.14) in some detail, using the LC+ approximation. With a starting state {p, f, c , c} m at time t 0 , we have (using eq. (5.20) in the LC+ version of eq. (5.14)) Here we have explicitly displayed the sum over parton indices in H I : The l and k dependent splitting operator is, from eq. (6.4), Here M c is the color matrix defined in eq. (6.5).
Now consider the no-splitting operator N LC+ (t, t 0 ). The color basis states are eigenstates of this operator, Using eq. (6.10), the integrand in the exponent can be written as an integral over momenta, a sum over flavors, and a sum over parton labels l and k: Here, χ(k, l, {c} m ) was defined in eq. (6.11) and N (k, l, {f } m+1 ) was defined in eq. (6.12). We see that the function λ({p, f, c} m , l, k, {p,f } m+1 ) that appears in the Sudakov exponent is almost the same as the integrand in the matrix element of H LC+ l,k (t). In fact  When we insert eq. (7.7) into eq. (7.1), we see that one can generate a splitting in standard Monte Carlo style by choosing the new momenta and flavors together with l and k with a probability proportional to λ({p, f, c} m , l, k, {p,f } m+1 ). This leaves a sum over the choices of colors, For each choice of {ĉ ,ĉ} m+1 , there is a color factor C and a statistical state vector {p,f ,ĉ ,ĉ} m+1 that is the input to the next splitting. In a computer program, one could imagine implementing the sum by summing the results returned by a splitting function that is called recursively. However, this is not really practical. Instead, one can perform the color sum Monte Carlo Style. One chooses {ĉ ,ĉ} m+1 with a probability ρ c , normalized to Then one assigns a weight factor to the splitting equal to Averaged over many trials, this reproduces the desired sums. The shower starts with a color weight factor of C 0 ({c , c} N 0 ) from the calculated color density matrix for producing an initial color configuration of N 0 final state partons at the start of the shower. At each splitting, we multiply by the weight factor w c from eq. (7.10). The shower evolves with successive splittings until it reaches some cutoff hardness value. At that point, let us say that we measure the expectation value of some color operator O c . Then this expectation value is We simply multiply the weight by this factor. In particular, making no measurement of color corresponds to multiplying by just the color overlap function {c } N {c} N . Thus if the final color is not measured, the complete generated event comes with a weight , where p c is the explicit power of 1/N c that is associated with the quantum amplitudes {c} N and {c } N as defined in section 5.11 and P (N ) is at least as big as the color suppression index I(N ) from eq. (5.31). Recall that the color suppression index either stays the same or else increases at each parton splitting.
One can use the behavior of the color overlap function to supply a hint about how to choose the probabilities ρ c . For numerical efficiency, one wants the dispersion of the weights not to be large. Consequently, it is sensible to make ρ c small for any choice of color configuration {ĉ ,ĉ} m+1 that increases the color suppression index. For instance, for the splitting shown in figure 7, the LC+ approximation retains the first and third terms on the right hand side; one would pick the color configuration shown in the first term with probability close to 1/2 and one would pick the color configurations shown in the third term with probability proportional to 1/N 2 c . We explore this idea further in appendix A using information from section 5.12 about how the color suppression index grows.
It is also possible to choose ρ c = 0 for any color choice {ĉ ,ĉ} m+1 that makes I(m + 1) greater than some predetermined limit I max on the color suppression index. This is an approximation beyond the LC+ approximation. We simply throw away configurations that have too many powers of 1/N 2 c . We explore this idea further in appendix B What is this operator? The leading color guess is that it is {c f } N {c f } N . However, this guess can be correct only in the leading color approximation because our color basis states are not exactly orthogonal to one another (and also the closed string states are not exactly normalized). Thus we need another set of basis states {c} N , ⊥ that equal our color basis states to leading order in 1/N c but are exactly orthonormal, The new states should be related to our basis states by a matrix equation where A is a unit matrix up to 1/N 2 c corrections. Eq. (8.1) implies that where There is a minimal solution to the equation A † A = G above. We let A be a real symmetric matrix equal to To find A exactly, one can diagonalize G and replace each eigenvalue λ by √ λ. Alternatively, the expansion allows one to write A as a power series in 1/N 2 c . We conclude that it is sensible to define Note that the numerator of w c involves matrix elements of A. It is presumably sufficient to calculate w c to leading order in 1/N 2 c . One possible string configuration is always One does not necessarily need to follow the LC+ shower all the way to a (1 GeV) 2 virtuality scale. One could use the LC+ shower for a few splitting steps, down to, say, a (100 GeV) 2 scale. Then one could assign a classical color string configuration {c f } N to the partonic state as outlined above. After that, one could use a leading color shower to get from the (100 GeV) 2 virtuality scale to the (1 GeV) 2 scale. This would be appropriate if the measurement to be made on the final state is only minimally sensitive to what happens at the softer scales.

Full color perturbatively
The shower evolution equation with full color has the form where U(t, t 0 ) obeys the evolution equation eq. (2.2), This differential equation can be solved iteratively in the form eq. (5.14) Here N LC+ (t 2 , t 1 ) is the no-splitting operator, It is well to recall here the essential point: the operator V LC+ (τ ) is diagonal in the standard color basis that we use, so that it is practical to calculate its exponential. Now, what if we want shower evolution with full color? Then we need This evolution equation is equivalent to which can be solved iteratively: (9.9) One can have any (small) maximum number of insertions of ∆H I (τ ) − ∆V(τ ) . For instance, to start with one might test whether one such insertion makes a significant difference. To make one insertion of ∆H I (τ 1 ) − ∆V(τ 1 ) , one would generate τ 1 at random. Then one would run an LC+ shower from the starting scale t 0 up to scale τ 1 . After that, application of ∆H I (τ 1 ) produces a weight and a new shower state, which is the starting point for an LC+ shower from τ 1 to the final shower time t f . Similarly, application of ∆V(τ 1 ) produces a weight and a new shower state, which is the starting point for a second LC+ shower from τ 1 to t f . The resulting values for the measurement function from the two second stage showers would then be summed. It is interesting to note the counting of logarithms in eq. (9.9). Suppose that the parton shower is used to calculate an observable O in which there is a large logarithm L, so that O has the form n c(n, 2n)α n s L 2n + n c(n, 2n − 1)α n s L 2n−1 + · · · . Suppose further that a shower with full color, generated by U(t f , t 0 ), correctly calculates all coefficients c(n, 2n) and c(n, 2n−1). Then the first term in eq. (9.9) correctly generates the coefficients c(n, 2n) since the LC+ approximation is exact with respect to color for the soft×collinear singularities. In the second term in eq. (9.9), the one insertion of ∆H I (τ 2 ) − ∆V(τ 2 ) generates a factor α s L by correcting the color content of a wide angle soft splitting. This factor multiplies factors α n−1 s L 2n−2 . Thus the second term makes contributions to the coefficients c(n, 2n − 1). The third and higher terms in eq. (9.9) contribute to coefficients c(n, j) with j ≤ 2n − 2 only. Thus to get the coefficients c(n, 2n − 1), one needs only the first two terms in eq. (9.9).
In general, by using one or two terms beyond the LC+ approximation in eq. (9.9), one can see whether the splitting operators beyond the LC+ approximation have an important influence on whatever observable is being investigated using the parton shower. One may expect that for many observables, the operators ∆H I (t) or ∆V(t) are not important. Including some factors of [∆H I (τ ) − ∆V(τ )] can test this hypothesis.

Soft gluon exchange phase
Up until now, we have presented the evolution equations for a parton shower in which the virtual splitting operator V is determined by the real splitting operator H I together with an assumption about the structure of V.
However, the structure assigned to V(t) leaves out an important physical effect: the wave function of two partons emerging from a hard interaction can accumulate an SU(3) phase factor U ≈ 1 + iφ(t)dt by exchanging a soft gluon. 10 Here φ is an operator on the quantum color state. In an inclusive cross section, this phase cancels between the virtual graphs for the bra state and for the ket state. That is, the phase cancels in a completely inclusive measurement because However, the color matrix φ(t) changes the color state, which influences the further evolution of the shower, so that the effects of the color phase do not cancel from all observables measured at the end of the shower. These effects can be important [11][12][13][14][15].
Recall the structure of eqs. (5.11) and (5.12) for V(t), which we can summarize in a shorthand notation as The phase is sometimes called the Coulomb gluon phase by analogy with the Coulomb phase in nonrelativistic quantum scattering. However, in a relativistic calculation in Coulomb gauge, only part of the phase comes from the Coulomb force between colored partons; the rest comes from the exchange of physically polarized gluons.
The probability conservation equation determines h through 11 This gave us h from H I (t), as in section 5.10. That is, we could find h without looking at any details of virtual graphs by simply knowing that the singularities of real emission graphs have to cancel with corresponding singularities of virtual exchange graphs. Previously, we made the assumption that φ = 0. However, when a soft gluon is exchanged between partons l and k, the graph has a part proportional to i times a hermitian matrix φ in the color space. This color phase is not zero. A simple calculation using the eikonal approximation, similar to other calculations in the literature (for example, ref. [11]), gives The phase depends on the relative velocities of partons k and l: Note that if either parton is massless, v kl = 1. There is a factor α s evaluated at a scale µ R that depends on the shower time and potentially depends on the momenta of partons k and l, depending the physical meaning of the shower time used in the parton shower (see appendix C). Finally, there is a color operator T k · T l . The operator T k inserts a color matrix T b on line k: if Ψ is the vector representing the color state {c} m before the virtual exchange as in eq. (2.4) or eq. (2.5), then the effect of the T k is to map Here b is the color index of the exchanged gluon. Similarly, T l inserts a color matrix T b on line l. Finally, we sum over the color index b, as indicated by the dot in T k · T l . The operator T k · T l maps {c} m into a linear combination of color basis states {ĉ} m . Let us denote the part of V(t) that comes from φ(t) by V φ (t). Imagine calculating the expectation value of some observable O by using a parton shower in which we calculate perturbatively in powers of V φ (t). We note immediately that contributions proportional to an odd number of factors of V φ (t) do not contribute to O because O is real and these contributions have an odd number of powers of i. However with an even number of powers of V φ (t) we have a factor ±(iπ) 2n = ±(−1) n π 2n and have a generally nonvanishing contribution. How big these contributions are depends on how sensitive the observable is to the color flow of the event. This is a question that deserves further study that is beyond the scope of the present paper. However, one can note immediately that the development of a parton shower depends crucially on the color structure of the shower state because this color structure determines the preferred emission direction for relatively soft gluons. Thus a sudden change in color structure caused by inserting two occurrences of V φ (t) at some early shower time t can have a substantial influence on the flow of momentum at the end of the shower. For instance, a gap in the rapidity of emitted partons could be created or could be filled in.
Can we say more about the color structure of the phase operator φ? We find that the color operator T k · T l applied to {c} m gives  12 The contribution to V φ from the first term in T k · T l in eq. (10.8) can be included in the LC+ approximation V LC+ for V. Then the remaining part of V φ becomes part of ∆V and can be treated perturbatively as in eq. (9.8).
The color phase operator φ depends on the parton masses. There are two places where the mass dependence might be important for the analysis of TeV scale processes at the LHC. First, in the later stages of showering, one can produce gluons with transverse momenta not too far above the b-quark mass. Sometimes such a gluon can split to b +b. Then the b andb are non-relativistic, so the mass dependence matters. Second, the hard process under investigation can produce top quarks or perhaps squarks, gluinos, and other very massive particles. In these cases, the particle masses matter.
Even though particle masses can matter, it is of interest to understand what happens when all parton masses can be neglected, so that v lk → 1. There can still be dependence on the parton labels k, l if the relation between the shower time and the renormalization scale in α s depends on the parton kinematics. Let us suppose that we neglect any such dependence. Then Since the whole shower is invariant under color rotations, we have used color basis states that are overall color singlets. Applied to color singlet states, we have k T k · T l = 0 (10.12) for any l. From this, we derive (10.14) In the first term in φ 0 , we sum over final state partons while in the next two terms we sum over the initial state partons. In either case, we have a color operator T i · T i , which is just the Casimir operator, with eigenvalue C A or C F depending on whether parton i is a quark or a gluon. Thus 15) where N F g is the number of gluons in the final state of {p, f } m while N I g is the number of gluons in the initial state and N F q is the number of quarks and antiquarks in the final state while N I q is the number of quarks and antiquarks in the initial state. When we apply −φ 0 to the bra state {c } m , we get exactly the opposite phase. Thus the term φ 0 in the color phase contributes nothing to the development of the shower and can be simply dropped. This leaves the single contribution φ ab , representing an effective double strength color exchange between the initial state partons. Notice that with the approximation used here there is no color phase for e + e − annihilation or for deeply inelastic scattering.
Using eq. (10.8) for the color operator T a ·T b , we obtain a leading color term that leaves the color state unchanged plus a term that changes the color state. The leading color term provides a phase that occurs but not the other, then there is a net phase that appears in every evolution interval until the color connection situation changes. These phase factors will tend to reduce the contribution of this sort of state to whatever observable is to be measured.
There are also contributions to T a ·T b that change the color configuration. These terms have the potential to change the energy flow in the final state by changing the color flow.

Conclusions
In general, a parton shower Monte Carlo event generator should generate contributions to the density operator in color space in which the color in the ket vector {c} m and the color in the bra vector {c } m are different. Standard parton shower generators based on evolution from hard splittings to soft splittings work in the leading color (LC) approximation, which is the leading order in an expansion in powers of 1/N 2 c . In the leading color approximation, only states with {c } m = {c} m occur.
In this paper, we have introduced the LC+ approximation, a generalization of the leading color approximation. Going beyond the leading color approximation inevitably involves sums over color choices. These sums can performed with Monte Carlo summation: selecting a choice at random according to a prescribed probability function. At the end of the shower each event comes with a weight.
The LC+ approximation has several nice features: • For each splitting, the leading soft×collinear singularity and the leading collinear singularity are treated exactly with respect to color. There is an approximation with respect to color, but it occurs only in wide angle soft splittings.
• Evolution can start with any state {c , c} m . Thus if one starts with a hard scattering process, one can fully use the color subamplitudes that multiply color states {c} m for the hard scattering, including in the calculation interference between ket states {c} m and bra states {c } m with different color configurations.
• The Sudakov factors are numbers. That is to say, the standard color basis states {c , c} m are eigenstates of the Sudakov operator. One does not have to exponentiate non-diagonal matrices in the color space.
• In fact, not only are the Sudakov factors numbers, but also there is a separate Sudakov factor for the ket state {c} m and for the bra state {c } m . This feature may prove useful for matrix element matching when the color structure of the amplitudes is treated exactly.
• With a simple extension of the formalism, one can include in the LC+ approximation the phase induced by exchange of soft gluons.
• The LC+ approximation is still approximate: remainder terms in the generators of shower evolution are left over. However, the remainder terms can be included in a perturbative calculation up to some order.
• The LC+ approximation can provide an efficient tool to sum large logarithms with full color at leading and next-to-leading log level for a certain class of observables if one uses the first perturbative correction as described in section 9.
The inclusion of weights in the shower generation has the potential to produce numerical problems. If the dispersion in the weights is large, then the number of events needed to calculate the expectation value of some observable with reasonable accuracy will be large. The total weight for an event is the product an initial weight, a final weight, and weights for individual splittings. By multiplying a large number of individual weights, one has the potential to produce large weights. For instance, (1+1/N 2 c ) N can be large if N is large. For this reason, it seems likely that there is a practical maximum to the number of splittings that can be generated using the LC+ approximation.
Fortunately, it is possible to turn the LC+ approximation off before it goes too far. In fact, we have found two ways to do that. First, one can simply run the LC+ algorithm for some number N max of splittings and then return to the LC approximation. For that, one must replace mixed states {c , c} m by diagonal states {c,c} m . We have presented a plausible model for doing that. Second, one can continue using the LC+ approximation, but not allow the generation of states with more than a certain amount of color suppression as measured by what we have called the color suppression index.
We are thus encouraged that the LC+ approximation can prove useful. The authors do not have at immediate hand a dipole based Monte Carlo event generator suitable for implementing the approximation. 13 However, the LC+ approximation is quite general and could, we think, be implemented in an existing generator.
Finally, one can ask if one really needs to go beyond the LC approximation. The answer is that we do not know if one needs to go beyond the LC approximation. The LC approximation is the first term in an expansion in powers of 1/N 2 c . We do not know if the contributions of order 1/N 2 c or beyond are important. One could imagine that these contributions are not important for some observables but are important for others. If that is the case, we do not know which observables are sensitive to 1/N 2 c effects. An implementation of the LC+ approximation would give us the means to investigate these questions.
In order to keep w tot from being large, it seems that the factors ρ c in the denominator of eq. (7.12) should be large, but there is a constraint from the normalization condition eq. (A.1): the probability budget should not be spent on configurations where it is not really needed. For that reason, one would set ρ c = 0 for color configurations {ĉ ,ĉ} m+1 for which C = 0.
One can also make ρ c small for color configurations that will, in the end, make where p c is the amplitude power defined in section 5.11 and P (m) is the color suppression power defined by this relation. Recall that P (m) ≥ I(m), where I(m) is the color suppression index defined in eq. (5.31). The color suppression index either stays the same or grows at each splitting step, so it is sensible to make ρ c small for any choice of color configuration {ĉ ,ĉ} m+1 that makes I grow.
We consider a g → g + g splitting for the case that k = l and that parton k is also a gluon. Look at M c , the numerator of C({ĉ ,ĉ} m+1 , {c , c} m , l m , k m , {f } m+1 ), as defined in eq. (6.5): Recall that C(l, m+1) selects color states in which partons l and m+1 are color connected. The operators t † l and t † k here insert the new gluon into the color states according to Here a † + (l) inserts gluon m + 1 just to the right of gluon l on its color string in {c} m and a † − (l) inserts gluon m + 1 just to the left. 14 This notation is explained in more detail in sections 7.2 and 7.3 of ref. [4]. The factor √ C F comes from the normalization of the color states, eqs. (2.4) and (2.5). Thus What does eq. (A.7) tell us? There are functions χ ± (k, l, {c} m ) and χ ± (k, l, {c} m ) that describe the color structure of the initial state. One of the functions χ ± (k, l, {c} m ) equals 1 if parton k is connected to parton l in the state {c} m . One of the functions χ ± (k, l, {c } m ) equals 1 if parton k is connected to parton l in the state {c } m . Potentially, parton k is color connected to parton l in both the bra and the ket state. If parton k is not color connected to parton l in either the bra and the ket state, then M c = 0. Then there is a color factor C F . Finally, there is a factor that describes the color structure of the state after the splitting. For instance, {ĉ ,ĉ} m+1 a † + (l) ⊗ a + (l) {c , c} m is 1 if gluon m + 1 is inserted to the right of gluon l in the ket state and in the bra state and this factor is zero otherwise. There are four possible ways to insert the new gluon, corresponding to the four terms in eq. (A.7). The first two terms correspond to parallel configurations while the last two terms correspond to crossed configurations. The fact that the crossed configurations come with minus signs tells us that we cannot use the coefficients in eq. (A.7) as Monte Carlo probabilities for the respective insertions.
We recommend defining the probabilities ρ c according to the status of gluon l in the color configuration {c , c} m . If gluon l is a healthy gluon, then the color suppression index remains unchanged if we chose one of the parallel configurations for the insertion of the new gluon, but if we chose one of the crossed configurations, then the color suppression index increases by 2. Thus we can chose parallel configurations with a probability close to 1 times an overall normalization factor and we can choose crossed configurations with a probability proportional to 1/N 2 c times the overall normalization factor. We define This is the counting factor that appears in the denominator of C in eq. (7.8).
If gluon l is a frail gluon, then the color suppression index remains unchanged for both parallel and crossed insertions of the new gluons. Therefore, we allow all of the color choices with equal probabilities:

B Avoiding too high color suppression
The LC+ approximation has the good feature that it can generate a shower starting from the color density operator terms {c} m {c } m with {c } m = {c} m . Such a state inevitably leads to {c } N = {c} N at the end of the shower and thus to a weight factor {c } N {c} N proportional to 1/N 2 c to some power. However, the LC+ approximation can generate color states with overlaps proportional to 1/N 2 c to a very large power. Suppose that we are content to limit the accuracy of shower evolution (within the LC+ approximation) to 1/N c to a certain power I max . For instance, one might set I max = 4. Then we can expand weight factors {c } m {c} m as a power series in 1/N c and not calculate contributions with too many powers of 1/N c . Furthermore, one can track the color suppression index I(m) of generated states. One can arrange not to generate any splitting that makes I(m+1) > I max .
Let us see how to avoid generating states with I(m + 1) > I max . Suppose that we have a color state with I(m) = I max . There may well be some frail gluons in the state. Within the LC+ approximation, these gluons can split. As we have seen in section 5.12, the splitting of the frail gluons leaves I(m) = I max and turns the frail gluons into healthy gluons.
There will likely also be some healthy gluons in the state. As we have seen, when a healthy gluon splits in the parallel configuration, the color suppression index is left at I(m + 1) = I max . However, when a healthy gluon splits in the crossed configuration, the color suppression index increases to I(m + 1) = I max + 2. This is beyond the limit that we want to maintain. Thus, when I(m) = I max , the splitting of a healthy gluon in the crossed configuration should not be allowed.
This analysis suggests that one should modify the probability function ρ c for a healthy gluon splitting when I(m) = I max by dropping the last two terms in eq. (A.8), so that one now uses This leads to a weight factor w c = 2C F /C A = 1 − 1/N 2 c . Since this procedure in any case drops color suppressed contributions once I(m) has reached I max , it makes sense to simply change the weight factor to w c = 1.
If one were to set the probability ρ c to generate a splitting of a healthy gluon in the crossed configuration to ρ c = 1, then one would very rarely generate an event with I(N ) > I max with a very large weight. We set = 0 so that no such events are actually generated. One should note, however, that then there is no control over contributions to the calculated observable at order greater than 1/N Imax

C Calculation of gluon exchange phase
To define the phase from soft gluon exchange, we consider a gluon with momentum q = (E, q) exchanged between possibly massive final state partons l and k. Using the eikonal approximation, and supplying an ultraviolet cutoff m 1 and an infrared cutoff m 2 , the amplitude for this process is θ(m 2 < | q| < m 1 )p l ·p k (p l · q + i )(p k · q − i )(q 2 + i ) . (C.1) Herep l andp k are the momenta of partons l and k as they enter the final state. We let It is convenient to defineθ by v l,x = cosθ . (C.18) (For massless partons,θ is half the angle between v l and v k .) Then the relation is 2| q| sin 2 (θ/2) + cosθ sin 2 (θ/2) ≈ Q 0 e −t . (C.19) We integrate over θ. Ifθ is of order 1, then the integration over θ is dominated by θ values of order 1. Ifθ 1, then the integration over θ is dominated by θ values of orderθ. Thus, after the θ integration, the effective relation between | q| and the shower time is 2| q| sin 2 (θ/2) ∼ Q 0 e −t . (C.20) Thus one possible relation between t and | q| is eq. (C.12), with C = 1/[2 sin 2 (θ/2)]. What about the renormalization scale, µ 2 R (t)? One might base this on the denominators of covariant Feynman diagrams instead of the energy denominators. Define D l (q) = (p l + q) 2 − m 2 l , wherep l is the final state momentum of parton l,p l = (E l , p s ) with E l = [ p 2 l + m 2 l ] −1/2 . The exchanged gluon is virtual, but let us consider the on-shell contribution with E = −| q|. Then D l (q) = 2p l · q = −2E l | q| + v l · q . (C.21) If we use −D l for µ 2 R we would take µ 2 R = 2E l Q 0 exp(−t). Similarly, we might take µ 2 R = 2E k Q 0 exp(−t). For the graph as a whole, might seem sensible as long as E l and E k are not too different.