The LPM Effect in sequential bremsstrahlung: $1/N_c^2$ corrections

An important question concerning in-medium high-energy parton showers in a quark-gluon plasma or other QCD medium is whether consecutive splittings of the partons in a given shower can be treated as quantum mechanically independent, or whether the formation times for two consecutive splittings instead have significant overlap. Various previous calculations of the effect of overlapping formation times have either (i) restricted attention to a soft bremsstrahlung limit, or else (ii) used the large-$N_c$ limit (where $N_c{=}3$ is the number of quark colors). In this paper, we make a first study of the accuracy of the large-$N_c$ limit used by those calculations of overlap effects that avoid a soft bremsstrahlung approximation. Specifically, we calculate the $1/N_c^2$ correction to previous $N_c{=}\infty$ results for overlap $g \to gg \to ggg$ of two consecutive gluon splittings $g \to gg$. At order $1/N_c^2$, there is interesting and non-trivial color dynamics that must be accounted for during the overlap of the formation times.


Introduction
Consider a high-energy gluon showering as it traverses a QCD medium, such as a quark-gluon plasma, via splitting processes such as gluon bremsstrahlung g → gg. At high energy, the formation time for a bremsstrahlung gluon becomes large and encompasses multiple scatterings with the medium, so that one must take into account the Landau-Pomeranchuk-Migdal (LPM) effect [1-3, 5-7, 9-11]. 1 It is possible for two consecutive splittings in the shower to have overlapping formations times. The corrections to an in-medium parton shower due to overlapping formation times are formally suppressed by a power of α s , and it has been of interest for many years to figure out exactly how significant such corrections are. 2 Such corrections were first analyzed in refs. [12][13][14] 3 for the case where one of the two overlapping splittings is relatively soft (with various other simplifying assumptions we will review later). Subsequently, the program of refs. [16][17][18][19][20][21] has been working toward analysis of the more general case where neither splitting is necessarily soft. However, the analysis of this more general case has used the large-N c approximation. The formalism for treating N c =3 is known in principle [22] 4 but may be challenging to implement numerically.
We'd like to know whether or not the N c =∞ overlapping formation-time calculations in the literature are a reasonable or poor approximation to the physical case of N c =3. In this paper, we investigate that question by calculating 1/N 2 c corrections to earlier N c =∞ results [16,18] for the effect of overlapping formation times on real double splitting g → gg → ggg. Our goal is to see whether those corrections, when extrapolated to N c =3, are large, small, or comparable to the purely parametric estimate O(1/N 2 c ) ∼ 10%. In this paper, other than going beyond the N c =∞ approximation, we will make the same sort of simplifying assumptions and approximations as in the earlier work of refs. [16][17][18][19][20][21]. We will assume that the medium is static and homogeneous on scales of the formation time and corresponding formation length. We also make the high-energy multiple scattering approximation and so take interactions with the medium to be described by theq approximation.
Before proceeding, we should clarify why the first corrections to N c =∞ are O(1/N 2 c ) instead of O(1/N c ). If one were to think inclusively about double splitting, then the g→gqq (pair production overlapping bremsstrahlung) rate would be an O(1/N c ) correction to the purely gluonic g→ggg rate because of the relative number of quark colors vs. gluon colors. However, though a calculation of g→gqq has not yet appeared in the literature without soft approximations, it may be computed using the same N c =∞ techniques that were used to compute g→ggg in refs. [16,18]. So computing g→gqq to leading order in the large-N c limit would give no information on the size of corrections to N c =∞ methods. Instead, we will focus in this paper exclusively on purely gluonic (overlapping) g→ggg. In theq approximation, the corrections to N c =∞ for the purely gluonic process will be O(1/N 2 c ). 5

Outline
In the next section, we first review the interesting, non-trivial color dynamics that take place for finite N c in calculations of overlapping formation time effects for g→ggg. We next discuss the N c →∞ limit and isolate the 1/N c and 1/N 2 c corrections to the effective "Hamiltonian" that describes medium-averaged evolution of high-energy gluons involved in the splitting process. In section 3, we start with what are called "sequential diagram" contributions to the g→ggg rate and show how to obtain analytic integral expressions for the 1/N 2 c corrections. The integrals can all be done analytically except for three time integrals, which later will be performed numerically. Section 4 fills in details about low-level N c =∞ formulas for applying theq approximation to different color combinations of the high-energy gluons involved in the g→ggg process. Section 5 generalizes the approach for sequential diagrams in section 3 to now also cover what are called "crossed diagrams." That completes the analytic work, and we move on to numerically evaluate the size of our 1/N 2 c corrections in section 6. We also discuss there the relation of our work to earlier work in a different context (relaxing the N c →∞ limit for single splitting rates g→gg that have not been integrated over transverse momentum). Finally, section 7 offers our conclusion.
We should clarify that, for simplicity, we will study 1/N 2 c corrections to only the subset of g→ggg processes that were studied for N c =∞ in refs. [16,18]. This leaves out, for example, direct g→ggg through a 4-gluon vertex, as opposed to a sequence of two 3-gluon vertices with overlapping formation times. Such direct 4-gluon processes have been studied in ref. [25] and found to be numerically small for N c =∞. Our study also leaves out effective 4-gluon vertices that appear in Light Cone Perturbation Theory from integrating out longitudinally polarized gluons in light-cone gauge. Their contribution even for N c =∞ has not yet been completed. (A calculation of their contribution in large-N f QED is included in ref. [20].) 2 Background: Color dynamics 2.1 Warm-up: The BDMPS-Z single splitting rate Throughout this paper, we draw diagrams for contributions to splitting rates using the conventions of ref. [16], which are adapted from Zakharov's description of splitting rates [9][10][11]. Fig. 1b gives an example for single-splitting (e.g. g → gg) in the medium. The high-energy particle lines shown in the figure are implicitly interacting with, and scattering from, the gluon fields of the medium, as depicted in fig. 2a, and the rate is implicitly averaged over the randomness of the medium. Such interactions with the medium change the color of each high-energy particle over time. At first, it may seem like calculating the rate would require a complicated analysis of the time dependence of the color of each such particle. value ofq. When making 1/Nc expansions in this paper, we treatq as fixed: we do not expandq in powers of 1/Nc. Our calculation of overlap effects for g→ggg in theq approximation therefore effectively involves only gluons. In standard discussions of large Nc for diagrams that involve only gluons, the expansion is an expansion in powers of 1/N 2 c [32]. In both cases, all lines implicitly interact with the medium. We need not follow particles after the emission has occurred in both the amplitude and conjugate amplitude because we will consider only p ⊥ -integrated rates. (See, for example, section 4.1 of ref. [16] for a more explicit argument, although applied there to a more complicated diagram.) Nor need we follow them before the first emission because we approximate the initial particle as on-shell. Only one of the two time orderings that contribute to the rate is shown above.
Fortunately, this is unnecessary for fig. 2a. 6 To get a flavor for the reason why, consider for a moment the extreme case where the medium itself is weakly-coupled. Then (to leading order in the coupling of the medium) the medium-averaged correlations of interactions with the medium are 2-point correlations, as shown in fig. 2b. Let's focus on one of these correlations, such as the green line connecting particles 1 and 3 in fig. 2c. Let T a n represent color generators T a (in the appropriate representation) that act on the color state of particle n. The interaction of particle 1 with the gluonic field of the medium comes with a factor of gT a 1 . The correlation of a pair of interactions of particles 1 and 3 with the medium then comes with a factor of (gT a 1 )(gT a 3 ) = g 2 T 1 · T 3 . But this operator is quite trivial because, by color conservation (after medium averaging), 7 the three high-energy particles in fig. 1b must form a color singlet, which 6 Here and throughout, we will only be considering rates which are fully integrated over the transverse momenta of the daughter gluons. Otherwise, the color dynamics is more complicated even for g→gg. See, for example, refs. [23,24]. 7 Without medium averaging, the color neutrality of the 3-particle state would not be conserved over time.
That's because the interactions in fig. 2a (via gluon exchange with the medium) may randomly change the color of just one of the three high-energy particles at a given moment, and exchanging one gluon with the medium turns a 3-particle color singlet into a 3-particle color octet. After medium averaging, however, the interactions with the medium must be correlated, such as in fig. 2b, and so color cannot flow out of the 3particle system since these correlations are instantaneous on the time scales relevant to splitting processes.
(In perturbative language, the medium-averaged correlator A a µ A b ν of background gluon gauge fields vanishes unless a = b.) The situation is analogous to translation invariance of a gas in thermal equilibrium: any particular configuration of the molecules is not translation invariant, but translation invariance is recovered after thermal averaging.  Figure 2. (a) g → gg but now depicting interactions with the medium. Here, each black line ending in a cross represents an interaction of a high-energy particle with the gluon field in the medium. (b) The medium average of those interactions in the case of a weakly-coupled medium. Here the black lines represent 2-point correlations of the medium interactions, which dominate for a weakly-coupled medium. (The 2-point correlations can be written in terms of correlations A a µ A b ν of the background gluon fields present in the medium.) The correlations are drawn as vertical in this time-ordered diagram because, in the high-energy limit, the correlation lengths in the medium are parametrically small compared to the length (time) scale of the high-energy splitting process. [Not shown but also present: short-time 2-point correlations between two medium interactions of the same high-energy particle.] (c) One correlation between particles 1 and 3 is highlighted.
where C i is the quadratic Casimir associated with the color representation of particle i. That means that T 1 · T 3 reduces to a simple fixed number in this context. (Specifically T 1 · T 3 = −C A /2 = −N c /2 in the case of g → gg.) Because of (2.1), we do not need to keep track of the dynamics of the individual colors of the three high-energy particles in order to calculate the rate for fig. 1. This conclusion can be generalized to strongly-coupled media as well when one describes medium interactions using theq approximation. See ref. [26] for the argument.
2.2 SU(3) color states for overlapping, double splitting Fig. 3 shows an example of a contribution to the rate for overlapping double splitting such as g → ggg. In the shaded region, the system has four high-energy particles (three in the amplitude and one in the conjugate amplitude). Again by color conservation, those four particles together must form a color singlet. Unfortunately, unlike the 3-particle case, color neutrality T 1 +T 2 +T 3 +T 4 = 0 is not enough to uniquely determine combinations like T i · T j which appear in correlations between high-energy particles' interactions with the medium. A similar uncertainty arises in more general arguments [26] in the context of theq approximation.
The source of this ambiguity is that there are many different ways one can make a color singlet out of four gluons (similar to how there are many ways to make a spin singlet out of four spin-1 particles in ordinary quantum mechanics). In SU(3), the color representations of two gluons can be combined as where the subscripts "s" and "a" indicate symmetric vs. anti-symmetric color combinations of the two gluons. We could make a color singlet out of four gluons by combining the first two gluons into any color representation R appearing on the right-hand side of (2.2), then combine the other two gluons into its complex conjugateR, and then combine the resulting R andR into a color singlet. This process is depicted schematically in fig. 4a, labeled "s-channel." These s-channel color states form a basis for all 4-gluon color singlet states. Alternatively, one may instead choose a "t-channel" or "u-channel" basis, as indicated in the figure. 9 In this paper, we find it convenient to work in the u-channel basis because of particle numbering conventions in earlier papers on overlapping formation times [16,18]. We will label the u-basis singlet states as |R u . Our initial basis for discussing 4-gluon singlets is then For the case where R = 8, we have to label whether each pair of the four particles formed the 8 by a symmetric (s) or anti-symmetric (a) combination, as distinguished in (2.2). As 9 For a variety of papers related to these constructions (and discussion of the color generalization of 6jsymbols to relate different channels), see, for example, refs. [22][23][24][27][28][29][30][31]. explained in the present context in ref. [22], 10 only a 5-dimensional subspace of (2.3) appears in calculations of overlapping formation times (e.g. fig. 3): 11 Soon, we will discuss how the color singlet state of the four gluons evolves in the subspace (2.4) as the gluons travel through the medium. But first, we wish to discuss the generalization from SU(3) to SU(N ).

SU(N ) color dynamics for overlapping, double splitting
For the sake of compactness, we will refer to the number of quark colors as N rather than N c in the rest of this paper. The generalization of the preceding discussion to N > 3 is that the tensor product (2.2) of two gluon colors becomes 12 A ⊗ A = 1 s ⊕ A a ⊕ A s ⊕ "10" a ⊕ "10" a ⊕ "27" s ⊕ "0" s , (2.6) where 1 is the singlet representation, A is the adjoint representation of SU(N ), and, for example, "27" means the SU(N ) representation that generalizes the 27-dimensional representation of SU(3). The scare quotes just mean that, though we quote the size of the representation for N =3, we really mean the corresponding representation of SU(N ). Note that there is one more term in (2.6) than in the original SU(3) product (2.2). This representation "0" of SU(N ) smoothly decouples and disappears as one approaches N → 3 from above. For SU(N ) with N > 3, there is a 6-dimensional (rather than 5-dimensional) subspace of color singlet states relevant to calculations of overlapping formation times, which is spanned by the basis [22,23] |1 u , |A aa u , |A ss u , |"10+10" u , |"27" u , |"0" u . (2.7) This generalizes (2.4). In this paper we will quote some results about 4-particle color singlet states from ref. [22], but we have found it convenient to use slightly different overall sign conventions for the definitions of the u-channel states (2.7). The details of the relation between our sign conventions here and those of ref. [22] may be found in appendix A.1.
In Zakharov's version of the BDMPS-Z calculation of single splitting rates, the problem is recast as two-dimensional quantum mechanics (in the transverse plane) with an imaginaryvalued "potential energy" V . Ref. [16] extended this picture, in the large-N limit, to calculations of overlap effects in double splitting, such as the contribution to the rate represented by fig. 3. The 4-gluon potential needed to treat the shaded region of fig. 3 for finite N was worked out in ref. [22] for theq approximation. The resulting 2-dimensional Hamiltonian for the 4-gluon evolution in the shaded region of fig. 3 was found to be 13 Above, symmetries have been used to reduce the 4-gluon quantum mechanics problem with transverse positions (b 1 , b 2 , b 3 , b 4 ) to an effective 2-particle quantum mechanics problem [16,18] written in terms of (C 41 , The P ij are the canonical momenta conjugate to the C ij , and E is the energy of the initial particle in the doublesplitting process. The x i represent the longitudinal momentum fractions of the four gluons. The underlined quantities in (2.8) represent 6 × 6 matrices (for N > 3) that act on the 6dimensional space of relevant 4-gluon color singlet states. The matrices S u and T u encode results for the action of T i · T j on this space in the u-channel basis (2.7), encoded as 14 13 See appendix A.1 for details of how the s-channel result of ref. [22] was translated to the u-channel version in (2.8).

N =∞ limit
In the N →∞ limit, (2.10) becomes (2.12) Unlike the case of finite N , the matrices S u and T u commute for N = ∞. It is therefore possible to find a new basis that simultaneously diagonalizes both matrices: 13) in terms of which the N =∞ limits (2.12) become (2.14) We will explain our naming convention for the basis states (2.13) shortly. We have dropped the subscript u on S N =∞ and T N =∞ just to keep our notation from becoming too cluttered.
Because S N =∞ and T N =∞ are both diagonal, the potential (2.8b), and so the Hamiltonian, does not mix the states (2.13) for N =∞. Each of these states propagates independently for N =∞, with non-matrix potentials given by using the corresponding eigenvalues from (2.14) in place of the matrices S u and T u in (2.8b). We will only encounter transitions between these color singlet states when we later investigate the O(1/N ) perturbations to S N =∞ and T N =∞ .
The motivation for the names |1 and |A ± in (2.13) should be clear enough. One may use the conversion matrices between bases given in appendix A.1 to see that the state |A × defined in terms of u-channel color singlet states is equivalent, in the limit N → ∞, to the combination |A aa s + |A ss s / √ 2 of s-channel basis states. Similarly, the state |1 × − is equivalent to the s-channel basis state |1 s , and |1 × + is equivalent to the t-channel basis state |1 t . So we may think of the cross "×" in the notation A × or 1 × ± as meaning that, for N = ∞, the state involves the representation R = A or R = 1 in a cross-channel different from our usual u-channel representation.
Later we will also use the definitions (2.13) of basis states when analyzing large but finite N . In that case the equivalences just discussed (and so the motivation for the notation) are not exactly correct. So, for N < ∞, one may also interpret the cross × in the colloquial sense of "crossed out": a warning that the motivation for the notation is no longer precise for those states.

An aside: Diagrammatic interpretation of basis states for N = ∞
We make a brief detour to present another way to characterize the basis (2.13) for N = ∞. This alternative characterization can offer insight and will be used for some detailed arguments in section 4.2, but is not strictly necessary for most of our calculation.
Refs. [16,18] discuss drawing time-ordered diagrams, such as fig. 3 and others on the surface of a cylinder, where time runs along the length of the cylinder. The large-N requirement that N =∞ diagrams be "planar" [32] can be translated to say that no lines should cross on the surface of the cylinder. So, for instance, fig. 3 can be drawn on the cylinder as in fig.  5, where we have numbered the lines during the 4-particle part of the evolution according to the convention of ref. [16], which for this diagram corresponds to identifying the longitudinal momentum fractions of the gluons as (x 1 , x 2 , x 3 , x 4 ) = (−1, y, 1−x−y, x). Correlation lines, such as the black lines drawn in fig. 2b (and also higher-point correlations), must also be part of the "planar" diagram and so must lie along the surface of the cylinder without crossing any other lines. As a result, for N =∞, there can only be correlations between high-energy particles that are neighbors of each other as one goes around the circumference of the cylinder. So, during the 4-gluon phase of the time evolution in fig. 5, the medium interactions of particle 1 can be correlated with those of particles 2 and 4 but not with particle 3. We will indicate this particular sequence as (1234). Any cyclic permutation, such as (2341), would be an equivalent designation, and so would the reverse order (4321) or its cyclic permutations. All that matters for discussing the interactions among the particles in large N is which of the four high-energy gluons are neighbors.
With this notation the color singlet states (2.13) may be identified as (see appendix A.2) when N = ∞. Above, the notation (ij)(kl) means that particles i and j are contracted into a color singlet and that particles k and l are also contracted into a color singlet.
In terms of the cylinder picture of fig. 5, representing states like (ij)(kl) requires two separate cylinders: one for each singlet pair. This is a useful convention because it corresponds naturally to the large-N topological principle that diagrams requiring handles are suppressed. Specifically, as a preview of what we will see later, fig. 6 shows one type of 1/N 2 correction to fig. 5. As time progresses during the 4-gluon part of the evolution, there is a 1/N suppressed transition from the (1234) color singlet state to the (12)(34) color singlet state, and then later another such transition to the (1243) color singlet state. In our notation (2.13), that's where each transition will be due to 1/N corrections to the Hamiltonian. Some examples of (2-point 15 ) correlations of medium interactions are shown by the black lines. In the language of large N diagrammatics, the resulting diagram (interpreted here to include the medium correlations shown) cannot be drawn as a planar diagram, which is why it is 1/N 2 suppressed. In general, there is a suppression by 1/N 2 for every handle needed to draw a diagram on a surface without crossing lines [32]. 16 2.6 1/N and 1/N 2 corrections to the potential We can now work out corrections to the N =∞ limit by expanding the original Hamiltonian (2.8) in powers of 1/N . The dependence on N appears only in the S u and T u matrices (2.10), which can be expanded in powers of 1/N . But we will want to express the result in the basis (2.13) of states that decouple in the N =∞ limit, not the original basis (2.7) used for presenting S u and T u . After that change of basis, matters is that no lines cross when the diagram and correlations are drawn on the surface. 16 See also Coleman's excellent "1/N " summer school lecture in ref. [33].
with S N =∞ and T N =∞ as in (2.14) and (2.17)

Sequential diagrams
The sample diagram we have been showing so far is called a crossed diagram [16] because two lines cross when it is drawn as in fig. 3 (as opposed to the drawing in fig. 5 of the same diagram on the cylinder). To study 1/N 2 corrections, it will be simpler to start with a different class of diagrams called sequential diagrams [18], shown in fig. 7. However, only the first diagram xyxȳ (and its complex conjugate and permutations) will generate 1/N 2 corrections. That's because, as discussed earlier, there is no interesting color dynamics for 3-particle propagation, which means that there are no finite-N corrections needed for those propagators provided one uses the value ofq appropriate for the desired value of N . (The same is true of 2-particle propagators.) Only the xyxȳ diagram in fig. 7 has a region of 4-particle evolution and so non-trivial color dynamics, denoted by the shaded region in fig. 8.

Set-up and allowed color singlet transitions
We now focus exclusively on the xyxȳ diagram. In fig. 8, our numbering of particles in the region of 4-particle evolution follows the same convention as refs. [16,18]. This diagram gives a contribution to the rate for overlapping double splitting g→ggg that is proportional to 17 .
Technically, integrating over all of the times (tx < ty < tx < tȳ) gives probability, not rate. We should integrate only over time differences, but that detail is unimportant for the present discussion.  . The above diagrams contributing to double splitting g→ggg are called the "sequential diagrams" in ref. [18]. As in refs. [16,18], the diagrams are individually named (xyxȳ, etc.) by the time order of the vertices. The relevant permutations referenced above are those permutations of the daughters x, y, and z ≡ 1−x−y that create distinct diagrams. Above (t x , t y , tx, tȳ) are the times of the four vertices in fig. 8 from left (earliest) to right (latest). The factors B y , t y |B x , t x and Bȳ, tȳ|Bx, tx represent the propagators for the 3particle evolution respectively before and after the shaded region of the figure. The factor Cx 41 , Cx 23 , tx|C y 41 , C y 23 , t y represents the propagator for the 4-particle evolution inside of the shaded region. There is a gradient ∇ (corresponding to a factor of transverse momentum) associated with each splitting vertex. We have not shown here other overall factors, including how those gradients are contracted together by helicity-dependent DGLAP splitting functions.
The only non-trivial corrections to N =∞ come from the color dynamics of the 4-particle propagator, which we now write as Our specification of the initial and final 4-particle states in the propagator (3.2) is incomplete: We will also need to specify what 4-particle color singlet states we start and end in. We find it convenient to rewrite (3.2) as is a 2-dimensional vector (with elements that are in turn 2-dimensional vectors in the transverse plane) encoding the transverse position state of the system at a given time; is the total duration of the 4-particle evolution; and λ y and λx label the initial and final 4-particle color singlet states for that evolution.
Ref. [22] explains that those initial and final singlet states are each |A aa u for the diagram of fig. 8. 18 A quick, graphical way to understand why is that (i) as far as color representations are concerned, everything to the left of the shaded region of fig. 8 looks like the u-channel diagram of fig. 4c with a gluon (R=A) for the internal line, corresponding to |A u ; (ii) 3-gluon vertices combine gluons anti-symmetrically via the group structure constants f abc , therefore specializing to |A aa u ; and (iii) there is no color dynamics for 3-particle evolution, which means that the interactions with the medium in the actual diagram of fig. 8 will not affect the correspondence with the color-contraction diagram of fig. 4c. A similar argument applies to everything to the right of the shaded region of fig. 8.
In terms of our N =∞ eigenstates (2.13), the initial and final color-singlet states of the 4-particle evolution are then So, we will be interested in 4-particle Green functions (3.3) where the initial state can be λ y = A + or A − and the final state can be λx = A + or A − . From the texture of the finite-N corrections (2.17) to the S u and T u matrices that appear in the Hamiltonian (2.8), we can now identify what 4-particle color-singlet transitions contribute to the 1/N 2 correction to the xyxȳ sequential diagram of fig. 8. As just discussed, the sequence of transitions must start and end with A ± . The transition sequences allowed by (2.17) are then Note that neither δS nor δ 2 T contribute to any allowed O(N −2 ) corrections for this diagram. There are no O(N −1 ) corrections to the diagram: neither δS nor δT produce a direct |A + → |A − or |A − → |A + transition. This is consistent with the fact that, for purely gluonic processes, corrections in a large-N analysis should appear in powers of 1/N 2 [32].
In passing, we note that the allowed transitions (3.7) can be written in the alternative language of (2.15) as Note that the last three sequences may be obtained from the first three sequences by exchanging (2↔3) particles 2 and 3. But the only thing differentiating particles 2 and 3 in the xyxȳ diagram of fig. 8 is their longitudinal momentum fractions y and z ≡ 1−x−y. This means that instead of calculating the contributions of all six sequences (3.8), one could, if desired, use only the first three sequences but then add (i) that result to (ii) the same calculation with the value of y changed to 1−x−y.

1/N perturbation theory for 4-particle propagator
Let G N =∞ λ represent the N =∞ 4-particle propagator for any of the N =∞ color singlet eigenstates of (2.13), indexed by λ. In perturbation theory in 1/N , then the transitions (3.7) correspond to O(N −2 ) corrections to the propagator of the form . (3.9) Above, t 0 = t y and t 3 = tx are the initial and final times of the 4-particle evolution (the shaded region) in fig. 8. The two O(N −1 ) perturbations to N =∞ evolution (caused by δT ) occur at intermediate times t 1 and t 2 , as depicted in fig. 9. Each λ ij designates an N =∞ eigenstate from (2.13). As discussed previously, the initial color singlet state λ 01 and the final color singlet state λ 23 must be |A + or |A − as in (3.6) and (3.7). δV (δT ) λ←λ represents the (λ, λ ) matrix element of the δT contribution to the potential (2.8b). The non-zero matrix elements are all the same because the non-zero matrix elements of δT in (2.17) are all the same: δV δV Figure 9. A depiction of the 2nd-order perturbative correction (3.9) in 1/N to 4-particle evolution.
The shading shows regions of 4-particle propagation where N =∞ propagators are used. The dashed lines represent insertions of the 1/N correction δV to the potential at intermediate times t 1 and t 2 , which are integrated over.
In order to focus on structure over details, and also to allow for later generalizations, we will find it useful to introduce some short-hand notation for (3.10) and also to distinguish the earlier-time and later-time insertions of δV in (3.9): (where we have used the fact that x 1 +x 2 +x 3 +x 4 = 0). Here, R (δT ) is a 2×2 matrix that mixes the two components (C 41 , C 23 ) of the vector ξ defined by (3.4). It does not do anything to the transverse position space in which each C lives except to contract the transverse indices, as in (3.10). If one wants to be explicit, one could think of the matrices R i shown in (3.11) as really being R i ⊗1, where the 2×2 identity matrix acts on transverse position space. However, in our discussion, we will not speak explicitly about the transverse space. So, for example, we will refer to ξ throughout this paper as a "2-dimensional" (rather than 4-dimensional) vector, and we will correspondingly refer to the matrices in (3.11) as the 2×2 matrices (3.12).
We will see that the integration over all intermediate transverse positions in our diagram can be performed analytically inq approximation. That will leave three time integrals (t 1 , t 2 , and ∆t ≡ t 3 − t 0 ) to later be performed numerically.
To continue, we need the structure of the N =∞ 4-particle propagators. In theq approximation, these are 2-dimensional harmonic oscillator propagators for a coupled set of two oscillators (C 41 , C 23 ). Adapting the notation of ref. [16,18], we will refer to the two complex normal-mode frequencies of this system as Ω (λ) ± and define the 2×2 diagonal matrix For our context here, we have introduced the subscript or superscript λ to indicate which color singlet-state (2.13) we are finding the N =∞ propagators for. Again adapting the notation of refs. [16,18], we will make a matrix a whose columns are the corresponding normal mode vectors: . (3.14) We will leave for later the details of exactly what Ω ± and a are for each N =∞ color singlet state λ. For now, we have enough to write out the structure of the harmonic-oscillator propagator, which is 19 17) and the prefactor 20

Integrating over ξ 1 and ξ 2
The integrals over ξ 1 and ξ 2 in the expression (3.9) for δ 2 G are related to Gaussian integrals and so may be done analytically. We find the results are more compact if we first combine the two integrals into a single Gaussian integral by defining a 4-dimensional vector 19 It is because we are working in the same basis (C41, C23) throughout the 4-particle evolution that the first and last terms in the exponent of (3.15) have the same matrix A λ . This is unlike the original N =∞ analysis of diagrams in ref. [16,18], where it was found more convenient to use a different basis at the two ends of the propagator. 20 For N =∞, calculations of individual time-ordered diagrams were ultraviolet (UV) divergent (even for tree-level processes), which was treated with dimensional regularization in ref. [19]. Those divergences, however, were associated with 4-particle evolution times ∆t → 0 and so with the vacuum limit of the 4-particle propagators G. For vacuum evolution, there is no interesting color dynamics, and it is color dynamics that our 1/N corrections describe. As a result, there will be no UV divergences in our calculations of corrections in this paper, which means that we do not need to use dimensional regularization and so may use the 2-transverse dimensional formula (3.18) for f λ .
from the two intermediate position vectors ξ 1 and ξ 2 . Similarly, define to be a 4-dimensional vector composed of the initial and final position vectors ξ 0 and ξ 3 for the 4-particle evolution. Then the expression (3.9) for δ 2 G, together with (3.11) for δV and (3.15) for G N =∞ , can be rewritten in the form where we define the 4×4 matrices Above, we use the shorthand notation The parameters j 1 and j 2 are dummy source term coefficients used to generate the two factors (3.11) of δV in (3.9) from the Gaussian integral appearing in (3.21). Doing that Gaussian integral gives 21 . (3.24)

Evaluating the xyxȳ diagram
We could now go through all the additional steps of (3.1) for evaluating the xyxȳ diagram, which involve taking gradients of the 4-particle propagator, including the initial and integrating analytically over the first and last vertex times t x and tȳ, and correctly keeping track of all the prefactors not shown explicitly in (3.1). Instead, we are going to use a trick to bypass all of that by realizing that we can adapt the final result of the same steps that were applied in the original N =∞ calculations of refs. [16,18]. The trick will be to cast the 4-particle propagator (3.24) for our 1/N 2 correction into the same schematic form as the 4-particle propagator originally used in N =∞ calculations. Let's first discuss the latter to introduce notation (X, Y, Z) that was used in refs. [16,18,19].
In the original N =∞ analysis of the xyxȳ diagram in ref. [18], there were two color routings that had to be considered, which in the language of our paper here correspond to taking the full 4-particle propagator Cx 41 , Cx 23 , tx|C y 41 , C y 23 , t y for this diagram in (3.1) to be either G N =∞ A − , corresponding to the two N =∞ eigenstates that appear in (3.6). The calculations in ref. [18] focused on the color routing called here A − = (1234), to which the result for the other color routing could be related by swapping the daughters y and z ≡ 1−x−y.
In evaluating the A − color routing, ref. [18] organized the calculation (following the method of ref. [16]) by writing the exponential piece of the corresponding harmonic oscillator propagator in the form 22 where the C ij -independent prefactor f is unimportant at the moment. The above equation just gives particular names to the entries of the matrices that in this paper we would call A A − and B A − : namely 23 where is a matrix that flips the vectors (C 23 , C 41 ) appearing in parts of (3.25) to the basis (C 41 , C 23 ) that we have used exclusively in this paper. For N =∞, particular formulas for the (X, Y, Z)'s were given in ref. [18], which also figured out how to write the final answer for the diagram in terms of the (X, Y, Z)'s. 22 Our (3.25) is not shown explicitly in ref. [18]. There the argument, in appendix E.2, proceeds by analogy with section 5.3 of ref. [16] and skips over this explicit formula. The analogous formula is eq. (5.41) of ref. [16]. 23 The relationship between (X , Y, Z) seq y and (X , Y, Z) seq x follows from eqs. (E.11-12) of ref. [18] and from our (3.33), which shows the relationship between our X here and the X in ref. [18]. Now compare the old N =∞ formula above to the contribution of a particular color singlet transition sequence λ 01 → λ 12 → λ 23 in (3.24) if we leave out the operation ∂ j 1 ∂ j 2 [· · · ] j 1 =j 2 =0 . The dependence on the C ij 's is then completely contained in the 4-vector z of (3.20) and so in the exponential factor of (3.24). Comparing this exponential factor with the one in (3.25), we see that it has the same form, except that the (X, Y, Z)'s for the N =∞ calculation are replaced by alternate versions, which we'll call (X,Ỹ ,Z), given by At the very end, we will also then need to restore the overall operation ∂ j 1 ∂ j 2 [· · · ] j 1 =j 2 =0 that we strategically ignored in order to relate the different calculations! Our starting point, the result of ref. [18] for the color routing A − = (1234), is 24 where Formulas for (ᾱ,β,γ), which represent various combinations of helicity-dependent DGLAP splitting functions, may be found in ref. [18]. The variables X seq y and X seq x are related to the variables X seq y and X seq x we introduced earlier in (3.25) by where the additional |M |Ω terms arise from the integration of the 3-particle propagators, as described in ref. [16]. Finally, the formulas for M i , Ω i , M seq f , and Ω seq f may be found in ref. [18]. These have to do with the 3-particle evolution (which has no interesting color dynamics), and they remain the same in our problem.
We now obtain the desired 1/N 2 correction to (3.30) by swapping the (X, Y, Z)'s to (X,Ỹ ,Z)'s and replacing the prefactor f (A − ) in (3.30) by the analogous non-exponential factors (and operations) in (3.24): . (3.34) We have summed over all color transition sequences in (3.7).
To forestall possible confusion, we should mention that the result (3.34) automatically includes the product u A aa |λ 23 λ 01 |A aa u = 1 2 (3. 35) of overlap factors of the initial and final 4-particle color singlet states (3.6) with λ 01 =A ± and λ 23 =A ± respectively. That's because the same set of factors, in the form of were already implicitly included in the N =∞ result (3.30) for the single color routing A − . 26

Correction to total sequential diagram rate
To get the 1/N 2 correction to the total sequential diagram rate, we need to (i) take 2 Re[· · · ] of (3.34) in order to include the correction to the conjugate diagramxȳxy, and (ii) add all 26 The language of color singlet state overlap factors does not appear in the original N =∞ calculation of ref.
[18]. But (3.36) is equivalent to the 1 2 in the factor 1 2 C 2 A discussed immediately after eq. (E.1) of ref. [18]. permutations of the three final gluons (x, y, z) which generate distinct diagrams. See fig. 10. Correspondingly, the total correction is The symbol "∆" on the left side of (3.37) is inessential to our present purpose and is included for the sake of consistency with the N =∞ discussion of ref. [18]. 27 ) Alternatively, one may use the discussion about y↔z after (3.8) to write where δ 2 A seq (x, y) is also defined by (3.38) except that the sum over allowed color sequences λ 01 →λ 12 →λ 23 in (3.34) is taken over only the first three sequences of (3.8). The appeal of the version (3.39) is just that it has a similar form to how N =∞ results have been previously presented [18]. 28

Color-representation dependent formulas
In order to use the preceding formulas, we need for each 4-particle color singlet state the corresponding normal mode frequencies and normal mode vectors for 4-particle evolution, with the vectors written in the (C 41 , C 23 ) basis that we have been using throughout. That is, we need formulas for the Ω (λ) ± and matrix a (λ) of eqs. (3.13) and (3.14). In this section, we 27 See, in particular, section 1.1 of ref. [18]. Because the 1/N 2 corrections to sequential diagrams come only from the xyxȳ diagram (and its conjugate and permutations), that distinction does not matter here. 28 See eq. (3.1) of ref. [18].
will present this information for all of our N =∞ eigenstates (1, A + , A − , A × , 1 × + , 1 × − ), not just the states that appeared in the xyxȳ transitions (3.7), because the other states will be useful later on in the evaluation of 1/N 2 contributions to crossed diagrams.
We will start from the results for the |A + and |1 color singlets. The others may be related to these using permutation symmetries, for which the alternate notation (2.15) for N =∞ color singlet states will be very useful.

|A − = (1234)
This is the canonical color state considered in the earlier, N =∞ papers such as [16,18]. A convenient summary of the relevant formulas for Ω ± and a can be found in eqs. (A. [21][22] and (A.27-30) of ref. [21], where our matrix a in the (C 41 , C 23 ) basis used here corresponds to the matrix called a y there. We note for later reference that these formulas all depend on the momentum fractions (x 1 , x 2 , x 3 , x 4 ) = (−1, y, 1−x−y, x) of the four gluons. So where Ω (λ) is the matrix defined in (3.13).

|1 = (41)(23)
In this note, the u-channel color singlet state |1 refers to the case where the particle pairs (41) and (23) are each contracted into a singlet. This yields simple normal modes in the (C 41 , C 23 ) basis. The 4-particle potential (2.8b) for N =∞ acts on the |1 state as The normal mode frequencies Ω ± and vectors (C 41 , C 23 ) ± are and Following refs. [16,18], the normal modes have been normalized so that where is the mass matrix whose inverse appears in the kinetic term of the Hamiltonian (2.8a) for the basis (C 41 , C 23 ) that we use here.
The corresponding modes (4.4) for (41) (23) were expressed in the (C 41 , C 23 ) basis. So, by making the same permutation 1 ↔ 3 to (4.4), we obtain normal modes for |1 × − in the (C 43 , C 21 ) basis: Since C ij = −C ji , we can convert to the (C 34 , C 12 ) basis (which we'll see is useful in just a moment) by negating (4.8) to get To convert to the (C 41 , C 23 ) basis used throughout this paper, now use the relation [16] 30 (4.10) to get 29 See the discussion of eqs. (5.16-18) of ref. [16]. Here we work in the basis (C41, C23) instead of (C34, C12), and so the indices 1234 there are relabeled 2341 here. 30 This relation comes from eq. (5.31) on ref. [16]. and Now permute the conversion (4.10) by 1↔4 and then use C ij = −C ji to get . (4.14) Applying this transformation to (4.13) then gives the normal modes in the desired basis:

|A + = (1324)
We can get this from the formulas for |A − = (1234) by similar permutation arguments.
Since C 32 = −C 23 , we may rewrite that as gives and Then use (4.14) to get

Crossed diagrams
We now turn to crossed diagrams for g → ggg. The canonical crossed diagram, to which all others can be related [16], is the xyȳx diagram shown in fig. 11.  Figure 11. The canonical "crossed" diagram. Particles in the (shaded) region of 4-particle evolution are numbered according to the convention of ref. [16].

Allowed Color Transitions
At the start of the shaded region of 4-particle evolution, the particles combine in the same way as for the sequential diagram of fig. 8, and so the initial 4-particle color singlet state is the same as before: However, the end of the shaded region is different: It is now gluons 1 and 2 that meet at a vertex. So, the final state is the s-channel version |A aa s rather than u-channel version (5.1). In this paper, we find it convenient to always stick to the definition (2.13) of our basis , which are defined in terms of u-channel singlet combinations. We need to figure out how to express our final (s-channel) color singlet state |A aa s in terms of this basis. The matrix that converts (for any N ) between the s-channel and u-channel versions of the original basis states (2.7) is given by [22,23,27] For our present purpose, the only piece of (5.2) that we need is Using (2.13) to convert to the basis states (1, A + , A − , A × , 1 × + , 1 × − ) that we use for our analysis in this paper, and then expanding in 1/N , For future reference, note that the overall sign of |A aa s is merely a phase convention choice for that state. Different choices of this sign convention must lead to compensating changes of sign in the rule for the diagrammatic vertex at the end of the 4-particle evolution in fig. 11. We will later discuss how to get the overall sign of our answer right without having to drill down into such details. 31 We may now using the initial and final singlet states (5.1) and (5.4), together with the textures of the perturbations δS, δT and δ 2 T of (2.17), to list all possible 4-particle color transition sequences that contribute to 1/N 2 corrections to the xyȳx diagram of fig. 11. They are listed in table 1.

2nd order in δV
We start by examining the first five lines of table 1, which are the cases that involve two insertions of perturbations δV in the 4-particle evolution. Schematically, these cases correspond to fig. 12, which is the crossed diagram analog of fig. 9. The formulas for these contributions to the crossed diagram xyȳx are basically the same as the formulas we found in section 3 for the sequential diagram xyxȳ except for some minor modifications. One modification is simply that we should use the color transitions given by the first five lines of table 1 instead of the sequential diagram transitions of (3.7). But there are other changes needed as well. 31 We did not have to think about the phase convention in our discussion of sequential diagrams because the initial and final color singlet states were both the same: |Aaa u. So changing sign convention |Aaa u → −|Aaa u would have no effect since the sign would appear twice in the calculation of the 4-particle evolution-once at the start and once at the end. Table 1. Allowed 4-particle color transitions at order 1/N 2 for the xyȳx diagram, along with (i) the associated δT , δS or δ 2 T factors, and (ii) the product of the initial and final color overlap factors λ i |A aa u and s A aa |λ f . Also shown is the product φ of (i) and (ii) relative to what it would be s A aa |A − A − |A aa u = 1 2 in the N =∞ calculation of the crossed diagram. The horizontal lines separate groups of processes that have to be handled differently: 2nd order in δV with two δT transitions; 2nd order in δV with a δT and δS transition; 1st order in δV with a T -based perturbation. There are no non-zero O(1/N 2 ) contributions at 0th order in δV . [Specifically, the 1/N 2 term in (5.4) for the final state |A aa s does not directly overlap the initial state |A aa u of (3.6).] δV δV Figure 12. The analog of fig. 9, now for the xyȳx crossed diagram.

Modification: (X, Y, Z)
The known N =∞ rate for the xyȳx diagram has a form similar to that quoted earlier for the A − =(1234) color routing of the xyxȳ diagram in (3.30). The xyȳx case is [16] 32 The (α, β, γ) are different combinations of helicity-dependent DGLAP splitting functions than those in the sequential case, and their formulas may be found in ref. [16]. The I n here have the same form as the I seq n of (3.31) except that the superscript "seq" should be removed from everything. However, the (X, Y, Z)'s are somewhat different from the (X seq , Y seq , Z seq )'s. In the original N =∞ calculation [16], they were defined so that the exponential factor in the 4-particle propagator was 33 Xȳ Yȳ Yȳ Zȳ (3.33). Explicit formulas for M i , Ω i , M f , Ω f may be found in ref. [16]. The pattern common to the presentations (3.25) and (5.6) of the sequential and crossed exponentials is that in each vector, the bottom C v ij is the one for which lines i and j come together at the corresponding vertex v of the diagram. It was the use of this convention that made the N =∞ rate formulas (3.30) and (5.5) for sequential and crossed diagrams have similar structure.
Similar to what happened for the sequential diagram, (5.6) for this crossed diagram just gives particular names to the entries of the matrices that in this paper we would call A A − and 32 Unlike N =∞ sequential diagrams, N =∞ crossed diagrams have only a single color routing. 33 See eq. (5.41) of ref. [16], with the caveat that, similar to our previous discussion of the sequential case, our Xy and Xȳ here do not contain the effects of the initial and final 3-particle evolution and are related to the Xy and Xȳ of ref. [16] by our eq. (5.7). B A − -the identically same matrices that were relevant to the case of sequential diagrams in (3.26). Here the relations are where S is the matrix from (4.10) that converts the (C 34 , C 12 ) basis into the (C 41 , C 23 ) basis: Like we did for the sequential xyxȳ diagram, we now want to put the exponential factor where S is again defined by (3.27). So, to compute (X,Ỹ ,Z)'s for the crossed diagrams, first compute A − B U −1 B as in section 3, then read out the (X,Ỹ ,Z) seq values using (3.29), and finally convert those values using (5.11) above.

Modification: the matrix R 2
The identification of the 2×2 matrices R 1 = R 2 = R (δT ) back in (3.12) was based on the fact that only δT transitions (3.7) were relevant for the sequential diagram xyxȳ. The same is true for the first three rows of table 1, which shows the allowed transition sequences for the crossed diagram xyȳx. We will refer to those first three rows as "δT δT " transition sequences.
In the next two rows of the table, however, the second transition of each sequence is instead a δS transition. We will refer to these rows as "δT δS" transition sequences. S u 34 If one removes all of the tildes, then the relations (5.11) also relate the N =∞ crossed and sequential formulas for (X, Y, Z), which can be verified from the formulas for (X, Y, Z) in refs. [16,18], once one uses (3.33) and (5.7) to isolate what we call the X 's from the X's. appears differently than T u in the potential (2.8b), and so its contribution to δV matrix elements will be different than those of the δT contribution (3.11). The non-zero matrix elements associated with δS are which can be written in the form of 1 2 ξ R (δS) ξ with So, the final rule is that we need to use , R (δT ) ) for δT δT transitions; (R (δT ) , R (δS) ) for δT δS transitions (5.14) in the construction (3.22a) of the 4×4 matrix U.

Final result for 2nd order in δV
With the preceding modifications, the final result for the first five rows of table 1 has the same relation to (5.5) that the sequential result (3.34) did to (3.30): In the context of the sequential xyxȳ diagram, we previously discussed that equality of initial and final 4-particle color singlet state overlap factors between (i) the calculation of  δV Figure 13. Like fig. 12 except with only one insertion of δV (or δ 2 V ) during the 4-particle time evolution.

A single δT or δ 2 T perturbation
We now turn to the last group of transitions in table 1, where there is only one δT or δ 2 T perturbation. Schematically, this corresponds to fig. 13. The analog of (3.9) is where t 0 = t y and t 2 = tȳ are the initial and final times of the 4-particle evolution. The analog of (3.21) is then Doing the Gaussian integral over ξ 1 yields the analog of (3.24): By comparison of the exponential in (5.21) to the N =∞ exponential in (5.6), and accounting for the change of basis (5.9), we find that for these processes Also, as before, the X 's are related to the X's by (5.7). The analog of (5.15) is then where the (X,Ỹ ,Z)'s are now those determined by (5.22) and the Φ λ 01 ,λ 12 is a normalization factor we discuss below. For simplicity, we will take R 1 = R (δT ) (5.24) in the definition (5.20) of U for all of the single δV processes summarized in table 2. Because we are using (5.24) for δ 2 T as well as δT perturbations, this leads to the first of two normalization issues. The operation ∂/∂j 1 [· · · ] j 1 =0 in (5.23) was constructed to introduce one factor of 1 2 ξ 1 R 1 ξ 1 (unexponentiated) into the calculation of the overall result. Based on (2.17), δ 2 T matrix elements relevant to the transitions in table 2 all have value 1/2N 2 . In contrast, the non-zero matrix elements of δT in (2.17) are all 1/ √ 2 N . To correct for this difference, our overall factor Φ in (5.23) will need to contain (among other things) a factor of for δT transition; Earlier, we explained that our starting point -the N =∞ result -implicitly contains an initial and final color singlet overlap factor (5.17), which equals the similar overlap factors s A aa |λ 23 λ 01 |A aa u needed for both δT δT and δT δS transition sequences. From the "color overlap" column of table 1, however, we see that these overlap factors are different for some of the other transition sequences. Our overall normalization in (5.23) will need to account for this difference as well. Putting that correction together with (5.25), the overall normalization correction we need in (5.23) is for δT transition;  Table 2. (They are the same as √ 2 N φ, where φ is the last column of table 1.) Note that the 1/N 2 behavior of (5.23) comes from two places: One factor of 1/N comes from the factor R 1 =R (δT ) (3.12) produced by the operation ∂/∂j 1 [· · · ] j 1 =0 , and the other comes from the values of Φ in table 2.

Correction to total crossed diagram rate
The 1/N 2 correction to the xyȳx diagram corresponds to the sum of the results of (5.15) and (5.23), each using the formulas for the (X,Ỹ ,Z)'s appropriate to that particular process [(5.11) or (5.22)] and each summed over the relevant entries of table 1. In order to connect with previous work, it will be convenient to give the name δ 2 C to the total ∆t integrand for the xyȳx diagram: To get the 1/N 2 correction to the total crossed diagram rate from eq. (5.27) for the xyȳx diagram, we need to follow the same steps as originally used for the N =∞ calculation in ref.  Figure 14. The sum of diagrams that define the quantity A(x, y) in ref. [16]. [16]. There, the total rate was organized by first summing over the diagrams represented by fig. 14 to get where we've used the same notation (A, B, C) as ref. [16]. 36 Finally, we need to sum over any remaining permutations of the daughters (x, y, z) that lead to distinct diagrams. Just as in the N =∞ analysis of ref. [16], this gives 37 6 Numerical results

Main results
To get results for 1/N 2 corrections, we numerically integrated over (t 1 , t 2 , ∆t) or (t 1 , ∆t). A short discussion of our numerical method is given in appendix B. For comparison, the N =∞ results from previous literature only require numerical integration over ∆t. 38 Dividing 1/N 2 corrections by the corresponding N =∞ result gives the relative size of the corrections. The relative size of 1/N 2 corrections to crossed and sequential diagrams for overlapping double splitting g→ggg are shown, respectively, in figs. 15 and 16 for N =3 (QCD). Our convention in these plots is to let y represent the energy fraction of the lowest-energy daughter, x represent the next lowest, and then z = 1−x−y represents the highest-energy daughter. So the plot region has been restricted by these conventions to y < x < 1−x−y. The 1/N 2 36 See eqs. (8.1-8.3) of ref. [16]. But, for the 1/N 2 term corrections being considered here, there are no additional "pole" terms, as previously discussed in footnote 20. For the same reason, it is also unnecessary to make the vacuum subtraction of eq. (8.4) of ref. [16]. 37 See eq. (8.1) of ref. [16]. 38 The N =∞ results for crossed and sequential diagrams were derived in refs. [16,18,19], but a convenient summary of results may be found in appendix A.2 of ref. [21]. (for N=3) Figure 15. The ratio, for crossed diagrams only, of (i) the 1/N 2 correction to (ii) the N =∞ result for the differential rate dΓ/dx dy for (the crossed diagram contribution to) overlapping double splitting g → ggg. We have used N =3 in this plot, but one may multiply the results by (3/N ) 2 to restore the N dependence of the 1/N 2 correction. Very tiny wiggles in the contour lines are an artifact of interpolation from a discrete set of numerical data points. We have left out y < 0.01 just to simplify the numerical effort that went into making this plot. The ratio goes to zero as y→0, as one may see from the later discussion of fig. 19 for a particular value of x.
corrections to sequential diagrams are very small: less than 1%. The corrections to crossed diagrams are substantially larger. The largest relative correction occurs at the apex of the triangular region, (x, y, z) = ( 1 3 , 1 3 , 1 3 ), where the correction is roughly 17%. We showed the separate crossed and sequential diagram ratios first because there is a subtlety to discussing relative corrections to the total rate (crossed plus sequential). 39 Fig. 17 shows a plot of the ratio total 1/N 2 correction total N =∞ rate (6.1) for g → ggg, but restricted to y > 0.1. Similar to fig. 15, there is a (local) maximum at the apex of the triangular region, where the 1/N 2 correction is roughly 17%. Unlike fig. 15, however, around y ∼ 0.1 the rate has started to grow with decreasing y. As we will explain, this small-y growth is an artifact of how we have so far chosen to look at the size of 1/N 2 corrections. Instead of showing the ratio (6.1), fig. 18 shows, for a particular value of x, the small-y behavior of (i) the N =∞ result for the total rate vs. (ii) the sum of the N =∞ result and the 1/N 2 correction. We've taken x = 0.37, which corresponds to the blue dashed line in fig. 17. 40 At small y, both the N =∞ results and the total 1/N 2 correction blow up 41 as 1/y 3/2 , and so we have followed the convention of ref. [18] and instead plotted π 2 xy 3/2 ∆ dΓ dx dy total (6.2) in fig. 18. The extremely tiny values of y plotted in the figure are unlikely to ever be relevant to any real-world physics because, at the very least, one needs yE T for our high-energy approximations. 42 Nonetheless, this figure is useful to understand the behavior of our formulas. The important feature is that the total N =∞ result crosses zero at y ∼ 0.01 (for this value of x). This is possible because ∆Γ/dx dy does not represent a rate; it represents the correction to a rate due to overlapping formation times (see section 1.1 of ref. [18] for explanation), 43 and a correction may be positive or negative. Where the N =∞ result vanishes, the relative size 40 There is nothing special about the specific choice x = 0.37. 41 For a hand-waving qualitative explanation, see section 1.4 of ref. [18]. (6.1) of the 1/N 2 correction to that N =∞ result will blow up to infinity, by definition. From fig. 18, however, one sees that there is very little difference between the N =∞ curve and the corrected curve for y < 0.1. In any case, in applications to energy loss and in-medium shower development, the 1/y 3/2 small-y behavior of overlapping double splitting g→ggg is canceled [21] by similar behavior of virtual corrections to single splitting g→gg, leaving behind doublelog divergences [12][13][14] that are independent of N . So the small-y behavior of fig. 18, and in particular the 1/N 2 corrections to the small-y behavior, are not of much physical interest. 44 Furthermore, the original motivation for the large-N approximation in this problem was as a tool to be able to study overlapping hard splittings (y ∼ x ∼ 1).
[18] explains why it is not meaningful to talk about such a "complete" rate of double splitting in an infinite medium. It has to do with the fact that one way to achieve g→ggg is via two independent single emissions g→gg that are abitrarily far separated in time. 44 The "double log" behavior referred to above arises from y −1 ln y behavior in the combined real and virtual rates, producing a double logarithm when integrated over y. This is in contrast to the more divergent y −3/2 infrared behavior shown in fig. 18 for the real rate by itself. The fact that the soft behavior is y −1 ln y when virtual corrections to single splitting are included has been known since the early work of refs. [12][13][14] on energy loss in the soft-y approximation. The lack of N dependence of those double-log results appears in their calculations as a special feature of the dynamics of the soft gluon emission limit. (An explicit calculation showing in detail the cancellation of y −3/2 divergences between real and virtual diagrams for the case N =∞ may be found in ref. [21], which is focused on generic-y results but also extracts their small-y behavior.) In principle, the best way to investigate the size of 1/N 2 corrections would be to quote the relative size of their effect on an (infrared-safe) characteristic of in-medium shower development. Since we do not have everything needed for that (such as 1/N 2 corrections for virtual diagrams), we interpret fig. 18 to mean that a reasonable proxy is the largest relative size of the 1/N 2 corrections to dΓ/dx dy for y values that are not small -namely, the roughly 17% correction at the apex of fig. 17.

More detail on small-y behavior of crossed vs. sequential
There are some interesting qualitative features about the small-y behavior of crossed vs. sequential diagrams. Fig. 19 shows various different elements that went into the previous small-y numerics of fig. 18.
First, as noted originally in ref. [18], the N =∞ crossed and sequential results individually behave like ln(y −1 )/y 3/2 even though their sum just behaves like 1/y 3/2 . Because of the log dependence of the individual contributions, we have chosen to plot the contributions to π 2 xy 3/2 ln(y −1 ) ∆ dΓ dx dy (6.3) here instead of the normalization (6.2) used for fig. 18. Notice that the absolute size of the 1/N 2 correction to sequential diagrams (green circles in fig. 19) is quite small compared to that for crossed diagrams (blue diamonds). This means that, in absolute terms, the total 1/N 2 correction is overwhelmingly dominated by crossed diagrams. [Sequential diagrams play a role in the size of relative corrections shown in fig. 17 because they affect the N =∞ denominator of (6.1) even though they do not noticeably affect the 1/N 2 numerator.] Fig. 20, where the normalization of the verical axis is returned to (6.2), shows that the crossed and sequential 1/N 2 corrections are different in another way as well. From this plot, we find numerically that the crossed diagram correction behaves like 1/y 3/2 for small y, with no ln(y −1 ) enhancement. This is why the corresponding blue diamonds in fig. 19 approach zero, due to the additional normalization factor 1/ ln(y −1 ) in that plot. In contrast, we find that the sequential diagram correction has an even milder dependence on small y, behaving like 1/y 1/2 . When integrated over y in applications, this means that the 1/N 2 correction from crossed diagrams will contribute to infrared (IR) divergences (similar to the IR divergences of the N =∞ results discussed in ref. [21]), but the 1/N 2 correction from sequential diagrams will be IR finite.  Figure 20. A plot for x = 0.37 of the y dependence of the different 1/N 2 contributions to ∆ dΓ/dx dy, multiplied by π 2 xy 3/2 as in fig. 18. Note that, in contrast to fig. 19, we have not divided by ln(1/y).

Comparison of size of 1/N 2 corrections to related work
In this paper, we have been focused on the problem of overlapping formation times for the double splitting process g→ggg. For simplicity, we have followed previous N =∞ work on this problem [16,18,21] and only considered rates ∆ dΓ/dx dy that have been integrated over the (small) transverse momenta p ⊥ of all three daughters. For technical reasons, 45 studying p ⊥ -integrated rates allows one to ignore what happens to any daughter after it has been emitted in both the amplitude and conjugate amplitude. It's the reason, for example, why the dynamics of the y gluon is no longer relevant after the first conjugate-amplitude (red) vertex in fig. 3 for xyȳx interference diagram. However, there is a different type of problem (not about overlapping formation times) where similar issues of 4-gluon color-singlet dynamics also arise: the un-integrated p ⊥ distribution dΓ/dx d 2 p ⊥ for single splitting g→gg in the medium. Unlike fig. 1 for the p ⊥ -integrated g→gg rate, one must instead follow the color dynamics for a time after the splitting has taken place in both amplitude and conjugate amplitude, corresponding to the shaded region of fig.  21. During this time, one must treat the color dynamics of the four gluons shown in the shaded region. 46 Refs. [23,24,34] have investigated how to treat this problem beyond the N =∞ limit. Unlike our work (on our different problem), their calculations approximate the 45 See, for example, the argument in section 4.1 of ref. [16]. 46 The color dynamics of the two daughters decouple after a time of order the formation time, often referred to in this context as the color decoherence time.  fig. 1a for the rate of single splitting g→gg, but here including later timeevolution of the daughters (shaded region) that must be included in order to study the p ⊥ distribution of the daughters. In the shaded region, the above interference term contains four gluon lines which, in the language we have used in this paper, requires treating 4-gluon color singlet dynamics in the medium.
trajectories of the high-energy particles as perfectly straight lines. So they only include color dynamics and not the dynamics of particle trajectories. In this rigid geometry approximation (also known as the "antenna" approximation), they are able to more easily treat finite or expanding media.
For our specific purposes here, ref. [34] is interesting because the authors explicitly calculate the 1/N 2 correction to the N =∞ limit (using a different approach to calculate 1/N 2 corrections than we have). In their numerics, they study a medium of length L with constant q. Their results will depend on L, and they study the case where the dimensionless ratio L q/E (which parametrically is the ratio of L to what the formation length L f would be in an infinite medium) is 0.55. In this particular case, they find that the 1/N 2 correction to the N =∞ distribution for g→gg can be as large as 16%. Direct comparison to our roughly 17% correction should not be taken too seriously, however, since (i) the process we study is very different, and (ii) their numerics hold quarkq F fixed as they vary N , whereas we hold gluonq A fixed.

Conclusion
With two caveats, we have found that 1/N 2 corrections to N =∞ results for overlapping double gluon splitting (g → ggg) can be as large as approximately 17% for N =3 (QCD). One caveat is the technical one explained in section 6.1 that measurements of relative corrections become meaningless when the leading (N =∞) answer goes through zero at small y ∼ 0.01, and so we have focused on the size of corrections for not-small y. Also, small-y emission is not the case where large-N techniques were necessary to simplify the problem, since previous work on overlapping formation times with a soft emission [12][13][14] (which included the effects of virtual emissions) was done without resorting to the large-N approximation. Our interest here has been to estimate the reliability of using N =∞ results specifically for the case where y is not small.
The other caveat is that, in this exploratory analysis, we have not included diagrammatic contributions to g→ggg that involve 4-gluon vertices nor, in Light-Cone Perturbation Theory, instantaneous longitudinal gluon exchange. Nonetheless, our provisional take-away is that the N =∞ limit taken in previous analysis is likely a moderately good approximation. Work on using the N =∞ approximation to answer the ultimate question about the effect of overlapping non-soft emission on in-medium shower development is ongoing (using results for N =∞ rates from ref. [21] and the framework suggested by ref. [17]).
Ultimately, a complete analysis of 1/N 2 effects on energy loss should also include calculation of virtual diagrams for g → gg, as discussed for N =∞ in ref. [21].
Calculating virtual diagrams through order 1/N 2 may also be interesting for better understanding soft radiative corrections to hard single splitting g→ gg. Such radiative corrections give rise to IR double logarithms [12][13][14] and sub-leading IR single logarithms. The single logarithms have been calculated for N =∞ (for infinite medium in theq approximation) in refs. [35,36]. We do not know for sure whether those single logarithms have any non-trivial dependence on N . It would be interesting to be able to explicitly check at order 1/N 2 .
Finally, to answer a question proposed in the introduction, we note that our roughly 17% corrections for N =3 are roughly consistent with (e.g. within a factor of 2 of) the naive, merely parametric guess of O(1/N 2 ) ∼ 10%.
However, (A.1) is not how u-channel states were defined in ref. [22]. There, t-channel states were first defined by 47 and then u-channel states were defined in terms of t-channel states by 48 where we will useū to denote the u-channel conventions of ref. [22]. Performing (A. The difference between our u-channel convention on the right-hand side of (A.1) and the convention on the right-hand side of (A.4) is The effect of this is to negate states that involves an anti-symmetric combination of particles 1 and 4. From (2.6), that's |A aa u and |"10 + 10" u out of our u-channel states (2.7). We can summarize this as The matrix that converts (for any N ) between s-channel and t-channel versions of the original basis states is given by refs. [22,23,27] as 49 47 See, for example, eqs. (2.14) vs. (2.15) of ref. [22]. 48 See eq. (2.9) of ref. [22]. 49 This specific conversion is adapted from table IV of ref. [22], which provides the entries of our (A.8) and whose last column provides the signs P in our (A.9). See footnote 13 of ref. [22] for discussion of how those results are related to refs. [23,27].
(A.8) (Note that V = V = V −1 .) As explained in ref. [22], the conversion between the s-channel and the u-channel basis of that paper is correspondingly Using (A.6), the conversion between the s-basis and u-basis in our paper here is given by (5.2a) with U = PVP, which equals (5.2b). For some purposes, it is useful to have the N =∞ limit of this matrix, which is (A.10)

A.2 Alternative descriptions for N =∞
Here, we will justify the N =∞ identifications made in (2.15). One relatively simple method is to start with a form where the potential V is written directly in terms of the transverse positions (b 1 , b 2 , b 3 , b 4 ) of the four individual particles instead of in terms of reduced variables such as the (C 41 , C 23 ) of (2.8b). This more direct expression is [22] 50 (This expression only assumes theq approximation and not that the medium is itself weakly coupled.) Now we can use the expressions (2.9) 51 for the T i · T j in terms of S u and T u . Then remember that for N =∞ our basis states (1, A + , A − , A × , 1 × + , 1 × − ) are simultaneous eigenstates of S u and T u with eigenvalues given by the corresponding diagonal entries of (2.14). All together, we can then write the explicit 4-particle potential for each of our basis states.
For example, for |A − , eq. (2.14) indicates that (S u , T u ) = ( 1 2 , − 1 4 ) for N =∞. Using (2.9), the corresponding potential (A.11) is then This is exactly the large-N behavior that we characterized as "(1234)," where each particle can interact only with its neighbors going around the cylinder. Doing the same thing with |A × , which has (S u , T u ) = (0, 0) for N =∞, gives which corresponds to what we called (1243). As another example, |1 × − has (S u , T u ) = (0, − 1 2 ) for N =∞, which gives where particles 1 and 2 interact only with each other, and similarly particles 3 and 4 interact only with each other. This corresponds to what we called (12)(34). One may similarly check all of the other identifications in (2.15).

B Numerical method
Here, we give a few comments about our numerical method. We have not tried to figure out the most efficient method but have adopted a brute-force approach implemented in Mathematica [37]. But there are issues, which we think may be useful to briefly summarize. First, because of round-off errors caused by subtractive cancellations, we find that we typically need to do intermediate calculations with much more than machine precision in 50 Specifically, see eq. (3.10) of ref. [22]. 51 See also footnote 14. order to succeed in a brute force approach. We therefore do our calculations with higher precision arithmetic in Mathematica.
Some of our formulas, like (3.34), involve derivatives such as ∂ j 1 ∂ j 2 [· · · ] j 1 =j 2 =0 . After some experimentation, we decided to implement these derivatives numerically as rather than doing the derivatives analytically (or using a more sophisticated numerical estimate of the derivative). 52 We let Mathematica handle the numerical evaluation of matrix inverses and determinants in our formulas.
We found that naive use of canned, adaptive integration routines took too much CPU time and caused a host of problems. Rather than working out how to tweak adaptive integration to do what we needed, we just did all of our integrals as simply as possible by using midpoint Riemann sums. For the ∆t integral (as opposed to the t 1 and t 2 integrals), we found it convenient to change variables as ∞ 0 d(∆t) f (∆t) = ∞ −∞ dz e z f (e z ). In practice, we then replace the infinite z integration region (−∞, +∞) by a finite region (z min , z max ) carefully chosen to cover everywhere the integrand is non-negligible [a choice which must be adjusted to study small values of y].
The cost of our brute-force method is that, because it is not adaptive, there are an annoying number of numerical approximations one must check to be sure that results are accurate, e.g. the size of in the numerical derivatives, the number of Riemann intervals for the (t 1 , t 2 , z) integrals, and the cut-offs (z min , z max ).
Finally, we should mention the strategies we used to attempt to avoid human error in our analysis and coding. The result (3.24) for the 1/N 2 correction δ 2 G to the 4-particle propagator was initially derived independently, and implemented numerically, by each of us in very different ways. One way was the method presented in the text. The other way did not use any tricks for packaging the transverse-position integrations into higher-dimensional vectors and matrices like eqs. (3.4), (3.19) and (3.22) but instead did the Gaussian integrals separately and explicitly, resulting in very long mathematical expressions for δ 2 G. Once the two methods agreed numerically, we then both switched to the method and code that most quickly produced results for δ 2 G, which was the method based on (3.24). Now using the same code for δ 2 G, we each independently produced or spot checked the various results in 52 If one wished to take the derivatives analytically, one could make use of small-j expansions such as and det(U −1 ) = det(U −1 0 ) 1 + j1 tr(U −1 0 R1) 1 + j2 tr(U −1 0 R2) + j1j2 tr(U −1 where we have promoted the 2×2 matrices R1 and R2 to 4×4 block-diagonal matrices by defining R1 ≡ R 1 0 and R2 ≡ 0 R 2 . However, this leads to more complicated formulas, which take extra CPU time to evaluate. this paper. Since we consulted with each other on general methods and development [e.g. equations like (3.29)], and helped each other to ferret out sources of numerical discrepancy, our work was not completely independent, but the most error-prone aspects of our numerical work were done independently.