The LPM effect in sequential bremsstrahlung 2: factorization

The splitting processes of bremsstrahlung and pair production in a medium are coherent over large distances in the very high energy limit, which leads to a suppression known as the Landau-Pomeranchuk-Migdal (LPM) effect. In this paper, we continue analysis of the case when the coherence lengths of two consecutive splitting processes overlap (which is important for understanding corrections to standard treatments of the LPM effect in QCD), avoiding soft-gluon approximations. In particular, this paper analyzes the subtle problem of how to precisely separate overlapping double splitting (e.g.\ overlapping double bremsstrahlung) from the case of consecutive, independent bremsstrahlung (which is the case that would be implemented in a Monte Carlo simulation based solely on single splitting rates). As an example of the method, we consider the rate of real double gluon bremsstrahlung from an initial gluon with various simplifying assumptions (thick media; $\hat q$ approximation; large $N_c$; and neglect for the moment of processes involving 4-gluon vertices) and explicitly compute the correction $\Delta\,d\Gamma/dx\,dy$ due to overlapping formation times.

When passing through matter, high energy particles lose energy by showering, via the splitting processes of hard bremsstrahlung and pair production. At very high energy, the quantum mechanical duration of each splitting process, known as the formation time, exceeds the mean free time for collisions with the medium, leading to a significant reduction in the splitting rate known as the Landau-Pomeranchuk-Migdal (LPM) effect [1,2]. 1 A longstanding problem in field theory has been to understand how to implement this effect in cases where the formation times of two consecutive splittings overlap.
Let x and y be the longitudinal momentum fractions of two consecutive bremsstrahlung gauge bosons. In the limit y x 1, the problem of overlapping formation times has been analyzed at leading logarithm order in refs. [4][5][6] in the context of energy loss of highmomentum partons traversing a QCD medium (such as a quark-gluon plasma). Two of us [7] subsequently developed and implemented field theory formalism needed for the more general case where x and y are arbitrary. However, we only computed a subset of the interference effects. Most significantly, we deferred the analysis of effects that require carefully disentangling (i) the computation of double bremsstrahlung with overlapping formation times from (ii) the naive approximation of double bremsstrahlung as two, consecutive, quantummechanically independent single-bremsstrahlung processes. In this paper, we compute the effects that require this careful disentanglement.
In the remainder of this introduction, we will first qualitatively discuss what effect overlapping formation times have on a simplified Monte Carlo picture of shower development, which will help us later set up the technical details of how to approach explicit calculations. We then give a more precise description of exactly what effects we calculate in this paper versus which are further deferred to later work. With those caveats, we present interim results for the example of real, double gluon bremsstrahlung g → ggg in the case of a thick medium. We wrap up the introduction with a very simple argument for the parametric size of the rate.
After the introduction, section II is given over to the calculation itself, where the most important issue will be the technical implementation of our method for isolating the corrections due to overlapping divergences. (Details of the calculation which do not involve new methods but instead closely follow those established in ref. [7] are relegated to appendices.) Final formulas for the case of a thick medium are summarized in section III in terms of a single integral, which may be computed numerically. Section IV offers our conclusion, including comments on the sign of the result.
A. Simplified Monte Carlo versus overlapping formation times

Overview
In this paper, we will ultimately present results by giving the correction ∆ dΓ/dx dy to double splitting, due to overlapping formation times, instead of giving the double-splitting rate dΓ/dx dy by itself. We begin by explaining why this is the physically sensible choice for the calculation that we do. shows the case where two consecutive splittings (colored magenta and green) happen to overlap. In all of these figures, we have exaggerated the transverse direction: for high-energy particles (and high-energy daughters), splittings will be very nearly collinear.
In order to simplify calculations in this paper, we are going to assume that the medium is thick-much wider than any of the formation times for splitting. Now consider an (approximately on-shell) high-energy particle that showers in the medium. That is, imagine that the medium is thick enough that there are several splittings in the medium, as shown in fig. 1a. In the situation pertaining to jet energy loss in quark-gluon plasmas formed in relativistic heavy ion collisions, this could apply to an energetic particle (not necessarily the primary energetic particle) that showers and stops in the medium. Imagine approximating the development of this shower by an idealized in-medium "Monte Carlo": Start with a calculation or model of the rate for each single splitting, assume that consecutive splittings are independent, and evolve through time, rolling dice after each small time increment to decide whether each high-energy particle splits then. Even for purely in-medium development, this description is simplistic and is not intended to describe the many effects that are included in actual Monte Carlos used for phenomenology. 2 We will refer to this idealized  FIG. 2: A pictorial version of the definition of ∆ dΓ/dx dy as the difference between (i) the doublesplitting rate (represented by the gray box) and (ii) the comparable rate given by idealized Monte Carlo (IMC) restricted to two splittings. Above, z ≡ 1−x−y.
calculation, based just on formulas for single-splitting probabilities, as the "idealized Monte Carlo (IMC)" result. The assumption that the splittings are quantum-mechanically independent is equivalent to saying that this idealized Monte Carlo treats the formation times as effectively zero. That is, the picture of fig. 1a is treated as fig. 1b.
We want to compute how to account for what happens when two of the splittings are close enough that formation times overlap, such as in fig. 1c, in which case the idealized Monte Carlo assumption that the splittings are independent breaks down. Let's ignore very-soft bremsstrahlung for now, since it is more-democratic splittings that naively dominate energy loss; the (nonetheless important) effects of very-soft bremsstrahlung have been treated elsewhere [4][5][6]. So imagine that in each splitting the daughters both carry a non-negligible fraction of the parent's energy. 3 In this case, the separation between splittings is on average large compared to the formation times, provided that α s (evaluated at the scale that characterizes high-energy splittings 4 ) is small. The chances that three or more consecutive (democratic) splittings happen to occur with overlapping formation times is then even smaller than the chance that two consecutive splittings overlap. So, to compute the effect of overlapping formation times, it is enough to focus on two consecutive collisions.
A cartoon of the corresponding correction is given in fig. 2. Let δH be the part of the Hamiltonian of the theory that includes the splitting vertices for high-energy particles. The first term on the right of fig. 2 represents a calculation of the double-bremsstrahlung rate dΓ/dx dy in medium, where we have formally expanded to second order in δH, even though the real-world situation may be that there are eventually many more splittings (such as in fig. 1c). The other terms, subtracted on the right-hand side of fig. 2, represent the result (dΓ/dx dy) IMC that one would obtain from an idealized Monte Carlo if one formally expanded that result to second order in the splitting probability. (More on what we mean by that in a moment.) Individually, both dΓ/dx dy and (dΓ/dx dy) IMC receive contributions from consecutive splittings that are separated in time by much more than the formation times, discussion) also attempt to heuristically model corrections to independent treatment of splittings [16,17]. In contrast to Monte Carlos used for detailed phenomenology, for some examples of theoretical insight gained from studying various characteristics of the idealized in-medium showers we focus on here, see refs. [11][12][13][14]. 3 As will be discussed later in section I D, the specific assumption here is (parametrically) that E daughter α 2 s E parent for each daughter. 4 For hard splitting in a thick medium, this is α s at scale Q ⊥ ∼ (qE) 1/4 . See, for example, the brief discussion in section I.E of ref. [7]. but those contributions cancel in the difference, as depicted pictorially in fig. 3. Indeed, individually, both dΓ/dx dy and (dΓ/dx dy) IMC as defined above would be formally infinite in the case of an infinite medium. (More on that in a moment as well.) In contrast, the result for is finite and depends only on separations the relevant formation time.

Some simple analogies
Since the above points are important, we will try to illuminate them with some analogies. First, as a warm-up, consider the decay of a particle in quantum mechanics. The generic way to compute the decay rate is to formally compute the probability P for decay, to first order in the process that causes the decay. One finds a result that formally grows proportional to the total time T as P = Γ T, from which we extract the rate Γ. But of course the probability of decay can't really be P = Γ T because that probability would exceed 1 for large enough T. Instead, the probability is P = 1 − e −ΓT . The formula P = Γ T is analogous to what we mean above when we say to formally compute or expand a result to a given order in perturbation theory. This example is a rough analogy to constructing an idealized Monte Carlo based on the single bremsstrahlung rate: Γ is analogous to the single bremsstrahlung rate, whereas the result P = 1 − e −ΓT is analogous to what you would get if you actually used a result for Γ in an idealized Monte Carlo (as opposed to discussing an idealized Monte Carlo result that had been "formally expanded" to some fixed order in perturbation theory). Let's now turn to an example that is more analogous to our double bremsstrahlung problem.
Consider the classical analogy of a very tiny device ("particle") that has a certain probability Γ 1 per unit time of emitting a flash of light. If we formally expand to first order in Γ 1 , then the probability of emitting exactly one flash of light in time T is (formally) P 1 = Γ 1 T. If we formally expand to second order in Γ 1 , then the probability of emitting exactly two flashes of light 5 in time T is P 2 = 1 2 (Γ 1 T) 2 . If we naively divided P 2 by T to get a rate, then we would awkwardly say that the rate for two flashes is Γ 2 = 1 2 Γ 2 1 T, which (unlike 5 Whether or not one treats the two flashes as distinguishable or indistinguishable is inessential to our analogy; and so the factors of 1 2 here and in the rest of the argument are inessential to the point we want to make.
FIG. 4: An example of a process in a corrected Monte Carlo that could easily be implemented if ∆ dΓ/dx dy is positive. Here, the 1 → 3 splitting, representing the inclusion of an additional possible process in the Monte Carlo with probability distribution ∆ dΓ/dx dy, would account for the correction due to the possibility of fig. 1c.
So the "rate for two flashes" (analogous to the dΓ/dx dy for double splitting) is not, by itself, a very meaningful quantity.
But now suppose that the device had the additional property that, for some interval δt after emitting one flash, the rate for emitting another flash was temporarily changed to Γ 1 + δΓ. We might then ask for the correction to the previous result. We could again formally expand to second order in flash rates (which are now correlated as just described) to find the probability P 2 in this new situation, which would roughly 6 have the form P 2 = 1 2 (Γ 1 T) 2 + Γ 1 δΓ δt T. If we divided by T to define a Γ 2 = 1 2 Γ 2 1 T + Γ 1 δΓ δt, we would again have something ill-defined as T → ∞. However, the correction to Γ 2 due to the change is perfectly well defined as ∆Γ 2 ≡ Γ 2 − 1 2 Γ 2 1 T = Γ 1 δΓ δt. In this analogy, Γ 2 is like our double-bremsstrahlung rate dΓ/dx dy; 1 2 Γ 2 1 T is like (dΓ/dx dy) IMC ; δt is like the formation time; and the correction ∆Γ 2 is like the ∆ dΓ/dx dy of (1.1).
To further illuminate the importance and relevance of the subtraction (1.1), we will present in a moment an analogy with the importance of similar subtractions in kinetic theory. But it will be useful to first discuss what one should do with (1.1) once one has calculated it.

Uses
What can one do with a calculation of the correction ∆ dΓ/dx dy? First note that its definition (1.1) is as a difference of two positive quantities. A priori, that difference might have either sign: negative if the effect of overlapping formation times suppresses the double bremsstrahlung rate and positive if it enhances it.
If the correction is positive, there is a relatively easy way to implement the correction to an idealized Monte Carlo simulation: Simply interpret ∆ dΓ/dx dy as the probability distribution of an additional type of local splitting process that "instantly" produces three daughters (instead of just two daughters) from one parent. This would allow for Monte Carlo showers such as fig. 4. In contrast, if ∆ dΓ/dx dy is negative, one has to work harder. (Examples in the context of relativistic heavy ion collisions: the Monte Carlo generators 6 We say "roughly" because, if one wants an exact answer, then there are boundary issues having to do with the end of time at t = T. These sorts of boundary issues will be important to address later on in our discussion but are not important for the purpose of this analogy.  5: A depiction of decay processes X → bb and X →bb for our first kinetic theory analogy, based on ref. [15]. JEWEL [16] and MARTINI [17] implement heuristic models for a reduction in multiplesplitting rates due to overlapping formation times. 7 ) In this paper, we focus on the calculation of ∆ dΓ/dx dy and will not pursue how to incorporate it into Monte Carlo. The earlier discussion of (uncorrected) idealized Monte Carlo was necessary, however, for the definition (1.1) of ∆ dΓ/dx dy. We will discuss the sign of the correction shortly, when presenting numerical results in section I C.

A kinetic theory analogy
The same issues that lead to focusing on ∆Γ rather than Γ (for higher-order processes) arise in kinetic theory problems as well. We start by briefly reviewing a kinetic theory example from the literature. Then we'll give an example more closely analogous to our double bremsstrahlung problem (and to our earlier flashing device analogy).
Consider a kinetic theory discussion of the following simple model, analyzed by Kolb and Wolfram [15] 8 in the context of baryogenesis in grand unified theories. The toy model has a stable, nearly-massless particle b and a massive, unstable boson X which can decay by both X → bb and X →bb, as depicted in fig. 5. If one writes a kinetic theory for these particles, one would include the processes of fig. 5 (and their inverses) in the collision term. Now consider additionally including the process bb → X →bb, depicted in fig. 6a. There is a problem of double counting because the Feynman diagram for bb → X →bb includes two different types of physical processes. One of these is the case where the intermediate X boson is approximately on shell, which we depict by fig. 6b. The other is what's left: the case where the X boson is off-shell, which we depict in fig. 6c with an asterisk on the label X. The first case ( fig. 6b) is already accounted for by solving kinetic theory using a collision term based of fig. 5, and so supplementing the collision term by the full result of fig. 6a would double count the case of on-shell X. Instead, the collision term should contain just (i) fig. 5 and (ii) the off-shell contribution of fig. 6c. How precisely does one define the "off-shell" part of the contribution? By rearranging the terms of fig. 6 to make it a definition, as in fig. 7. This subtraction is (crudely) analogous to our problem's subtraction (1.1). The fact that fig. 7 and not fig. 6a is the correct thing to use in the collision term in Kolb and Wolfram's problem is analogous to our discussion of adding ∆ dΓ/dx dy rather than dΓ/dx dy to a Monte Carlo description originally based on single-splitting rates, as in 7 On a related note: See ref. [18] for a discussion of implementing negative-weight corrections in the context of vacuum Monte Carlo. 8 See in particular the discussion surrounding eq. (2.3.12) of ref. [15] FIG. 6: A depiction of (a) the process bb → X →bb separated into contributions (b) where the X is on-shell plus (c) the rest. The wavy double-line cut in the middle of (b) is used to indicate that fig. 6 to show the subtracted quantity that should be used in the collision term of a kinetic theory analysis. fig. 4.
We now give a kinetic theory example that is somewhat more analogous to the current problem. Ignore the LPM effect, but consider a kinetic theory description of a QED plasma that includes the leading-order, 2→3 process for bremsstrahlung, depicted by fig. 8. What if we now want to systematically include higher-order processes in the collision term? Consider in particular the 3→5 process depicted in fig. 9a. Just like fig. 6, this process contains two types of contributions. Including the contribution with an on-shell intermediate line, depicted by fig. 9b, would be double counting. Instead, one should only add to the collision term the remaining piece, fig. 9c, which is defined as the difference between figs. 9a and 9b. Here, figs. 9a, b, and c are analogous to our problem's dΓ/dxdy, [dΓ/dxdy] IMC , and ∆ dΓ/dxdy respectively.
These are not perfect analogies (they do not involve the LPM effect), but we hope they help illuminate the importance of the subtraction (1.1).

B. What we compute (and what we do not)
The preceding work [7] developed most of the formalism we will need for carrying out calculations and then (in approximations reviewed below) computed the subset of contributions to the double-bremsstrahlung rate depicted in fig. 10. It is very convenient to alternatively represent these contributions as in fig. 11, where the upper (blue) part of the diagrams depict a contribution to the amplitude and the lower (red) part depict a contribution to the conjugate amplitude. Ref. [7] referred to these as the "crossed" contributions to the rate because the interior lines are crossed in the representations of fig. 11. In both figs. 10 and 11, we explicitly show only the high-energy particles; the (many) interactions of those high-energy particles with the medium are implicit. (See the introduction of ref. [7] for more

description.)
In this paper, we will now evaluate the diagrams of figs. 12 and 13, which we refer to as "sequential" contributions because the two bremsstrahlung emissions happen in the same order in both the amplitude and conjugate amplitude. To compute the desired correction ∆ dΓ/dx dy to double bremsstrahlung due to overlapping formation times, we will need to subtract from our results the naive calculation of double bremsstrahlung as two consecutive, quantum-mechanically independent splitting processes, as in (1.1). In the last two diagrams (xxyȳ and xxȳy) of figs. 12 and 13, the x and y bremsstrahlung processes do not overlap in time. We will see later that these diagrams roughly, but not quite, match up with the idealized Monte Carlo calculation. Figuring out how to correctly compute the difference will be the main new technical development required for this paper.
As discussed in the preceding work [7], it is possible to set up the formalism in a quite general way that would require both highly non-trivial numerics and a non-trivial treatment of color dynamics to implement, but one can proceed much further analytically by making a few additional approximations. Though the methods we will discuss in this paper can be applied more generally, it behooves us to keep things simple in this first treatment. So we will follow ref. [7] when it comes to explicit calculations, by making the following approximations.
• We will assume that the medium is static, uniform and infinite (which in physical terms means approximately uniform over the formation time and corresponding formation length).
• We take the large-N c limit of QCD to simplify the color dynamics.  • We make the multiple-scattering approximation to interactions with the medium, appropriate for very high energies and also known as the harmonic oscillator orq approximation. 9 Also as in ref. [7], we will focus on the case where the initial high-energy particle is a gluon (and so, in large-N c , the final high-energy particles are also all gluons), as the resulting finalstate permutation symmetries make for fewer diagrams to consider so far. However, there is a downside. In the case of gluons, one must also consider the 4-gluon interaction, which gives rise to additional interference contributions. Examples are given in fig. 14. Because the calculation of these additional diagrams would distract from the main point of this paper, which is how to treat the sequential diagrams in the calculation of ∆dΓ/dx dy, we will leave fig. 14 for future work [21]. It turns out that these contributions are small whenever at least one of the three final gluons is soft, so the still-incomplete results that we derive here will have some range of applicability. But we will need fig. 14 for a complete calculation for the case of arbitrary x and y.
Another problem that we defer for another time is the change in the single-bremsstrahlung rate due to virtual corrections, such as the one shown in fig. 15. This has been worked out in the limiting case y x 1 in the context of leading parton average energy loss in refs. [4][5][6] and is related to anomalous scaling of the effective medium parameterq with energy.
Finally, we should mention that the relative transverse momentum of daughters immediately following a double-splitting has been integrated in our results. Though we will make some qualitative comments regarding transverse momenta later on, we have not calculated transverse momentum distributions. (If we did not integrate the double-splitting rate over transverse momenta, the calculation would be much harder, including the necessity of accounting for decaying color coherence effects [22] occurring after the last splitting time in, for example, each of the diagrams of fig. 10. See appendix F for a brief discussion.)

C. Preview of Results
Numerical results are given in fig. 16 for the sum of the crossed and sequential contributions to the correction ∆ Γ/dx dy for real double gluon bremsstrahlung from an initial gluon. 10 As mentioned above, the effects of 4-gluon vertices (such as in fig. 14) have not yet been included. We do not expect these to be important when one of the final gluons is soft, and so we can already draw the conclusion from fig. 16 that sometimes the correction ∆ Γ/dx dy is positive and sometimes it is negative, with the corresponding implications for the relative ease or difficulty of Monte Carlo implementation discussed earlier at the end of section I A. However, as a pragmatic matter, note from the figure that one final gluon has to have about 10 times smaller energy than the other two in order to get a negative correction. So, if one did implement a correction to real double bremsstrahlung in Monte Carlo, the most important effects for shower development (that are not simply absorbable into the running ofq from soft emissions as discussed in [4][5][6]) may correspond to the cases of positive correction, which is the more straightforward case to implement. However, the most interesting region is where all three final gluons carry substantial momentum fractions, and for that case we will need those 4-gluon vertex contributions, which we have left for later work.
In the special limiting case of y ≪ x 1, the numerical results approach This result is negative, as in the corresponding region of fig. 16. In a moment, we will give a simple general argument about why ∆ dΓ/dx dy could be expected to scale parametrically 10 As noted in the figure caption, the three final state gluons are identical particles. Here and throughout the paper, our ∆ Γ/dx dy is normalized so that rates are (formally) related to differential rates by We say "(formally)" because the correction ∆Γ to the total rate based on ∆ dΓ/dx dy is infrared divergent. like 1/xy 3/2 , as above. This scaling suggests that a smoother way to plot our final results, even outside of the y x 1 limit, is to pull out a factor of 1/π 2 xy 3/2 from the answer. Fig. 17 presents our results that way. One needs to see a larger range of y to verify the (somewhat slow) approach to the limit (1.2), which is shown in fig. 18 for a particular, small fixed value of x.  fig. 16 but removing a factor of 1/(π 2 xy 3/2 ) from the answer. The apex (x, y) = ( 1 3 , 1 3 ) of the triangular region above corresponds to π 2 xy 3/2 ∆ dΓ/dx dy = 1.05 [in units of We should note that the soft limit (1.2) for the correction ∆ dΓ/dx dy to the double real gluon bremsstrahlung rate cannot be extracted from early work [4][5][6] on soft multiple bremsstrahlung because that work used approximations that are only valid for the sum of (i) real double gluon emission with (ii) virtual corrections to single gluon emission and are not valid for (i) and (ii) individually. 11 Also, our result (1.2) depends on diagrams that were 11 Though refs. [4][5][6] give final results that are integrated over y, it is possible to extract the integrand and so extract something one might be tempted to call dΓ/dx dy. In particular, it's possible to identify the parts of the calculation that correspond to diagrams representing real double emission from those representing not evaluated in refs. [4][5][6], such as the xȳxy diagram of figs. 10 and 11.
The integral formula we will derive for ∆ dΓ/dx dy in the general case is a complicated expression that is painstaking to implement. In Appendix A, we provide, as an alternative, a relatively simple analytic formula that has been fitted to approximate fig. 17 very well.

D. Why 1/xy 3/2 ?
Before moving on to the details of calculations, we give one crude but simple explanation of why the rate for overlapping emissions has parametric dependence (dΓ/dx dy) overlap ∼ 1/xy 3/2 when y < x < z ≡ 1−x−y.
First, let's review some features of the LPM effect for the usual case of single bremsstrahlung, with x representing the less energetic of the two daughters. In the QCD virtual corrections to single emission. See Appendix F of ref. [7] for a particular example showing how one of our diagrammatic contributions to dΓ/dx dy exactly matches, in the y x 1 limit, a corresponding piece of the calculation of Wu [6]. However, this cannot be done for all contributions to dΓ/dx dy because the methods of earlier work such as [5,6] assumed that the colors of the highest-momentum particle in the amplitude and the highest-momentum particle in the conjugate amplitude jointly combine into the adjoint representation (A ⊗ A → A). This is a valid assumption (at leading-log order in the y x 1 limit) only if one adds together real and virtual emission diagrams. See Appendix F of ref. [7]. case, the formation time for this bremsstrahlung process is (for thick media) parametrically 12 Each formation time offers one opportunity for splitting, with probability of order α s . So the rate for emission of a daughter with momentum fraction of order x is Now consider double bremsstrahlung with y < x. From (1.3), the y formation time will be shorter than the x formation time. The probability that an emission with momentum fraction of order y happens during the x formation time t form,x is then given by 12 See, for example, Baier [23] for the QCD formation time expressed in terms ofq, or the discussion leading up to eq. (4.15) in the review [24] (where µ 2 /λ ∼q and ω = xE). Here's one quick argument: The least energetic daughter is the most easily deflected. If it picks up transverse momentum of order Q ⊥ from the medium during the splitting process, the angular separation between it and the other daughter will be of order the ratio of its transverse and longitudinal momenta: θ ∼ Q ⊥ /xE. Over time t, the transverse separation of the x daughter from the trajectory the parent would have followed then grows to be of order b ∼ θt ∼ Q ⊥ t/xE. The quantum coherence of the emission process ceases when this separation can first be resolved by the transverse momentum Q ⊥ of the splitting process-that is, when b ∼ 1/Q ⊥ . This condition defines the formation time t form . By definition ofq, the average Q 2 ⊥ ∼qt. Combining the above relations using t = t form yields (1.3). 13 A quick note on infrared cut-offs: The power-law divergence of (1.4) as x → 0 may at first sight look like the LPM effect is causing an enhancement of the splitting rate in this limit, but the LPM effect is always a suppression. A useful way to understand this is to rewrite (1.4) as where τ el is the mean free time between elastic collisions with the medium and N ∼ t form /τ el is the number of elastic collisions during the formation time. The first factor α s /τ el is (parametrically) the rate in the absence of the LPM effect, in which case each collision offers an independent opportunity for bremsstrahlung. The 1/N factor is the LPM suppression, and the above analysis assumed N 1 (i.e. t form τ el ) so that we could, for instance, take Q ⊥ ∼ √q t form in the preceding footnote. When this assumption fails completely (t form τ el ), the rate is simply the unsuppressed rate Γ x ∼ α s /τ el .
Multiplying (1.5) and (1.7) gives the rate for overlapping x and y emissions: If the presence of the y emission has a significant effect on the x emission (or vice versa), then the correction ∆ dΓ/dx dy defined by (1.1) would then also be This is indeed the parametric behavior (1.2) of our result in this paper. It turns out that, separately, the crossed and sequential diagrams are each logarithmically enhanced, for y x 1, but the logarithmically-enhanced contributions cancel each other in the total (1.9). Appendix B discusses how the logarithmic enhancement of various individual contributions can be related to collinear logarithms from DGLAP evolution and fragmentation, and how there is a Gunion-Bertsch-like [25] cancellation of those logarithms in the total result. This is not to say that collinear logarithms are never relevant, just that they do not appear at the order α 2 s /xy 3/2 of (1.9). They do appear, notably, in calculations [4-6] of energy loss. In that case, the effect of the leading contribution (1.9) to energy loss from double bremsstrahlung is canceled by the leading contribution from virtual corrections to the single bremsstrahlung rate, and so it is the sub-leading behavior of each which becomes important.
Finally, note using (1.3) that the total probability (1.6) of emitting a y during the formation time of x is of order So one would need to deal with resumming multiple O(y) emissions during the x formation time if interested in y α 2 s x. We will not treat that case in this paper. 14 (When discussing the mathematical behavior of our results, we will nonetheless discuss some absurdly small values of y, such as in fig. 18. Our results for such small y directly apply only to the formal limit that the α s associated with splitting is arbitrarily small. 15 FIG. 19: The xxyȳ interference, showing longitudinal momentum fractions x i and our notation for the vertex times.

II. THE CALCULATION
In this section, we turn to the specifics of calculating the sequential diagrams of fig.  13. The first interference diagram in the figure, xyxȳ, can be evaluated by mostlystraightforward application of the methods of the preceding work [7]. We will leave the details for later, in section II B and appendix E. The only new subtlety of this diagram, as opposed to the crossed diagrams evaluated in ref. [7], is that there are two inequivalent ways that color can be routed in the large-N c approximation, which must be accounted for appropriately.
To begin, we will instead focus on the other explicit diagrams of fig. 13, xxyȳ + xxȳy, as these are the diagrams that involve the most significant new issue: careful attention to the subtraction of (dΓ/dx dy) IMC from dΓ/dx dy.

A. 2 Re(xxyȳ + xxȳy) vs. idealized Monte Carlo
Start by considering xxyȳ. Our convention for labeling time in this diagram is shown in fig. 19. Following the philosophy of the preceding paper [7], we may interpret the evolution between vertex times to be described by 2-dimensional non-relativistic non-Hermitian quantum mechanics, with the appropriate number of particles. In the case of fig. 19, it is a (i) 3-particle problem for t x < t < tx, (ii) 2-particle problem for tx < t < t y , and (iii) 3-particle problem for t y < t < tȳ. Also following ref. [7], we may use symmetry to reduce the N particle problems to (N −2) particle problems, leaving us with (i) 1-particle, (ii) 0-particle, and (iii) 1-particle problems respectively.
At this point, we could implement the methods of the preceding paper [7] to turn fig. 19 into explicit equations and then start turning the crank. We do this in appendix C for the sake of concreteness. However, many of these details can be sidestepped if we take a slightly looser approach, which is how we will proceed here in the main text.
The important point is that, for the xxyȳ diagram of fig. 19 (and similarly for the xxȳy diagram in fig. 13), the time interval (t x , tx) of the x emission does not overlap with the time interval (t y , tȳ) of the y emission. So the times of these events are unrelated-there is no reason the y emission cannot occur a very long time after the x emission. In the formalism, this is because evolution during the intermediate time interval tx < t < t y corresponds to a problem with effectively zero particles and so does not have any time dependence: the xxyȳ interference contribution does not care how far apart tx and t y are. This is consistent with interpreting this diagram as representing two consecutive splittings that are completely independent from each other, as in an idealized Monte Carlo calculation.

A crude correspondence
In particular, the xxyȳ + xxȳy cases (plus conjugates) in figs. 12 and 13 are related to the idealized Monte Carlo contribution of fig. 20, where a first splitting, E → (1−x)E and xE, is followed later by an independent second splitting, (1−x)E → zE and yE. In what follows, it will be convenient to introduce the longitudinal momentum fraction y of the y daughter with respect to its immediate parent in the second splitting, That is, the y daughter has energy yE = y(1 − x)E. In the language that we use for labeling diagrams, single-splitting rates are given by 2 Re(xx) = xx +xx, with xx depicted in fig.  21. The correspondence here between double-splitting diagrams and Monte Carlo is, crudely speaking, where T formally represents the amount of time between the first splitting and the end of eternity (a regularization that we will soon have to treat more carefully As in the earlier discussion in section I A, the infinite quantity T will only make a temporary appearance along the way towards finding the correction ∆ dΓ/dx dy to idealized Monte Carlo. Let's now use some more explicit formulas for the single-splitting rates appearing on the right-hand side of (2.3). In notation similar to that of the preceding paper [7], 16 corresponding to twice the real part of the xx diagram for single splitting, shown in fig. 21. (2.5) Here, P (x) is the relevant (helicity-averaged) DGLAP splitting function. B is the single transverse position associated with the "effectively 1-particle" description of this process. The subscripts at the end of Bx, ∆t|B x , 0 E,x are a reminder (which will be useful in a moment) of what energy and branching fraction should be used in the calculation of the propagator for B. The integration variable ∆t in (2.4) represents the time difference tx − t x in fig. 21. We will later give more explicit formulas for (2.4) in the multiple scattering (q) approximation, but for the moment there is no reason not to be general. To simplify notation in what follows, we will rewrite (2.4) generically as (2.7) Putting the integral form (2.6) of the single bremsstrahlung rate into the right-hand side of (2.3), we can write the idealized Monte Carlo result for double bremsstrahlung as As we verify explicitly in appendix C, the above idealized Monte Carlo contribution corresponds to the result for the diagrams xxyȳ + xxȳy +xxȳy +xxyȳ = 2 Re(xxyȳ + xxȳy), with one critically important difference. To explain that difference, look at the piece of (2.8) corresponding just to xx followed by yȳ [as opposed to the full 2 Re(xx) = xx +xx followed by 2 Re(yȳ) = yȳ +ȳy]: The actual result from the xxyȳ diagram is the same except for the constraints on the time integration in fig. 19: 19).
(2.11) (There are only three time integrals on the right-hand side because of our focus on computing rates rather than probabilities in this paper, and our associated assumption of time translation invariance over relevant time scales.) If we are sloppy, we might be tempted to conclude that these two types of time integrations are the same by (i) rewriting the right-hand side of (2.11) as (2.12) (ii) noticing that the integrand in (2.10) does not depend on t y − tx, just as previously discussed concerning the propagation of the effectively 0-particle intermediate state in the diagram for xxyȳ; and so (iii) replacing the d(t y − tx) in (2.12) by T to get the left-hand side of (2.11). The last step is not quite correct, however, because the maximum allowed duration of the intermediate time interval (tx, t y ) in fig. 19 will depend, for example, on how much of the available time T is being used up by the following time interval (t y , tȳ), since tȳ must also occur before the end of eternity. The integrand in (2.10) will fall off exponentially when the duration tȳ − t y of the last interval becomes large compared to the y emission formation time, and so that interval's effect on the sloppy result is only an edge effect, changing the right-hand side by subtracting an amount of order T 0 . Such an edge effect would be negligible (relatively speaking) as T → ∞ were we interested only individually in the formal result for xxyȳ or the corresponding idealized Monte Carlo contribution (results which are individually physically nonsensical as T → ∞). However, we are instead interested in the difference ∆ dΓ/dx dy, for which the leading O(T) pieces cancel, leaving behind the O(T 0 ) corrections, which then have a sensible T → ∞ limit. We could attempt to proceed by figuring out how to somehow consistently carry through the calculation with a sharp cut-off on time, as we were sloppily considering above. This introduces many headaches. 18 The path of least confusion is to introduce a smooth and physically-realizable cut-off. One regularization method is to imagine a situation where the strengthq of medium interactions very slowly falls off at large times to reach the vacuum. For example, takeq (t) =q 0 e −t 2 /T 2 . (2.14) The prescription that then arises in the T → ∞ limit turns out to be reasonably intuitive and straightforward to implement.

The prescription
Here, we will state the prescription and use it to find ∆ dΓ/dx dy. We leave to appendix D a more thorough justification, based on slowly-varying choices ofq(t) like (2.14).
Let τ x and τ y be the (instantaneous) times of the x and y emission in the idealized Monte Carlo treatment of fig. 20. When comparing to the the calculation of 2 Re(xxyȳ + xxȳy), the prescription is to identify these Monte Carlo times with the midpoints of the corresponding time intervals (t x , tx) and (t y , tȳ) of the emissions. That is, Now remember that the T in the idealized Monte Carlo formula (2.2) was just an attempt to regularize the integral Instead of giving that integral a name T, let's just step back and write this integral in place of T. We can then combine it with the ∆t integrals in the more explicit idealized Monte Carlo result (2.10) to recast (2.10) as dΓ dx dy IMC corresponding 18 One headache is that there is usually difficulty with radiation fields whenever you allow a charged particle to suddenly appear or disappear. Another is that we've just seen that edge effects will contribute to our answer, but if a single splitting takes place right at the edge of time, the presence of the edge will affect its rate. For instance, if t y in fig. 19 occurs a third of a formation time before the end of time, then there is no room for tȳ − t y to stretch out as far as one formation time. So you would not reproduce the single splitting rate at the edge, making comparison with idealized Monte Carlo problematic.
with the ∆t's and τ 's defined here by (2.9) and (2.15). In contrast, the result of the xxyȳ diagram is given using the integral on the right-hand side of (2.11), The difference, which contributes to ∆ dΓ/dx dy, corresponds to tx<tx<ty<tȳ dtx dt y dtȳ · · · − τx<τy dtx dt y dtȳ · · · (2.19) which may be reorganized as (2.20) The integrand in (2.17), however, depends only on the ∆t's and not separately on τ y , so the dτ y integration in (2.20) is trivial and leaves (2.21) Letting the notation [∆ dΓ/dx dy] xxyȳ represent the difference between xxyȳ and the corresponding piece of the idealized Monte Carlo calculation, (2.17) becomes Treating xxȳy similarly (as well as the conjugatesxxȳy +xxyȳ), The important thing about these integrals is that there are no large-time issues because ∆ dΓ/dx dy cares only about the cases where the formation times overlap. The integrals above are infrared (∆t → ∞) convergent. The ultraviolet (∆t → 0) is another issue, which we will address later, similar to the small-∆t divergences encountered for the crossed diagrams in the preceding paper [7]. Now specialize to the multiple scattering approximation, where the quantum mechanics problems are harmonic oscillator problems. As reviewed in our notation in the preceding paper [7], 19 the single-splitting result is then specifically (2.24) For simplicity of notation, we specialize now to the specific case treated in this paper, where all of the high-energy particles are gluons. In that case, the complex harmonic oscillator frequency appearing above is The earlier result (2.23) is then where The integrals in the result (2.26) have logarithmic divergences associated with ∆t → 0: We will see in section II B that this divergence cancels a similar small-time divergence of the first diagram in fig. 13, xyxȳ (plus its conjugate).
As discussed in the preceding paper [7], one must be careful about pole contributions from ∆t = 0 even when 1/∆t divergences cancel. We will defer these pole contributions to section II C.

Permutations
So far, we have only talked about subtracting away one contribution to the idealized Monte Carlo result-the one shown in fig. 20. What about the others? A complete summary of the ways that Monte Carlo can produce three daughters (x, y, z) is shown in fig. 22. Let's focus on the case where all high-energy particles are gluons. Then all the idealized Monte Carlo contributions in fig. 22 are related by permutations of (x, y, z), and they will be subtracted in the corresponding permutations of (2.26). The first Monte Carlo contribution in fig. 22 is subtracted from 2 Re(xxyȳ + xxȳy), as we have been discussing. The second is subtracted from the x ↔ z permutation 2 Re(zzyȳ + zzȳy). The third is subtracted from the x ↔ y permutation 2 Re(yȳxx + yȳxx).
Note that the xxyȳ diagram of fig. 19 maps into itself under y ↔ z, given the identity of our final state particles. So we should not include this permutation when we sum up all the diagrams. Later, though, it will be slightly convenient if we arrange our notation to pretend that y ↔ z is a distinct permutation. So, looking ahead, define Φ(x, y, z) by but also as Φ(x, y, z) + Φ(y, z, x) + Φ(z, x, y) + Φ(y, x, z) + Φ(x, z, y) + Φ(z, y, x). (2.31)

B. Discussion of xyxȳ
We now turn to the remaining sequential diagram: the first diagram (xyxȳ) of figs. 12 and 13. Evaluating this diagram is mostly an exercise in applying the methods of the previous paper [7], and we leave most of the details to appendix E. However, there is one new issue that we touch on here in the main text: the different ways one may route color in the xyxȳ diagram in the large-N c limit.

Color routings
Like our discussion of xxyȳ above, the xyxȳ diagram of fig. 13 technically remains the same if one permutes y ↔ z. However, in this paper we are working in the large-N c limit, shows the corresponding color flow for an example of medium background field correlations (black) that gives a planar diagram (and so leading-order in 1/N c ). In our notation, this interference contribution could be referred to as either xyxȳ 1 or xzxz 2 .
and, in that limit, there are two distinct color routings of the xyxȳ diagram which are not individually y ↔ z symmetric. We show these two large-N c color routings in figs. 23 and 24, which we will refer to as xyxȳ 1 and xyxȳ 2 respectively. In the figures, we follow the convention of the preceding paper [7] of drawing our large-N c , time-ordered diagrams on a cylinder. In large N c , correlations of high-energy particles' interactions with the plasma only exist between high-energy particles that are neighbors as one goes around the cylinder, which is why in this context the diagrams of figs. 23a and 24a represent different diagrams. Note that xyxȳ 1 and xyxȳ 2 are related by y ↔ z, and so we could also call them xzxz 2 and xzxz 1 respectively. The distinguishing difference between the calculations of the two color routings will, in the language of the preceding paper [7], be the assignments of the longitudinal momentum fractions x i for the 4-particle part of the evolution t y < t < tx. Going around the cylinder, the first routing xyxȳ 1 has (as labeled in fig. 23) whereas the second routing xyxȳ 2 has (as labeled in fig. 24), Because this last assignment is identical to the one used for the canonical diagram analyzed in ref. [7], we will focus on the evaluation of xyxȳ 2 to simplify comparison with previous work. Then, if we define the full result that we want from xyxȳ plus its distinct permutations, including all distinct color routings, is We mention in passing that the reason we could sidestep discussion of color routings for the xxyȳ diagram analyzed in section II A is because in that case there was no interval of 4-particle evolution and so no distinction like (2.32) vs. (2.33). Unlike 4-particle evolution, the choice of ordering of the x i for 3-particle evolution makes no difference to the calculation since 3 high-energy particles are all neighbors when drawn on our time-ordered cylinder.

Result and small-∆t behavior
In notation similar to that used in the preceding paper [7], the final result for xyxȳ 2 is where formulas for the various symbols are given in appendix E and ∆t represents the intermediate time interval tx − t y . This formula is identical to that for xyȳx in ref. [7] except for the addition of a superscript "seq" on some symbols (standing for "sequential" interference diagram as opposed to crossed), the bars on (ᾱ,β,γ), and subscript labelsȳ there replaced byx here. These modifications are explained in the appendix.
The small ∆t limit of the integrand is also discussed in the appendix, with result dΓ dx dy xyxȳ 2 = α 2 s P (x) P (y) As in ref. [7], the 1/(∆t) 2 divergence may be eliminated by subtracting out the vacuum calculation (which must total to zero when summed over all interference processes). If we now add in the conjugate diagramxȳxy 2 , that leaves This result is invariant under y ↔ z and so will be the same for the other color routing xyxȳ 1 . In total, then For real ∆t, this indeed cancels the small-∆t behavior of (2.28), as promised. So our final ∆t integrals for the sum of all sequential diagrams ( fig. 13) will be convergent.

C. Pole Contributions
As mentioned earlier and as discussed in the previous paper [7], one needs to be careful about the cancellation of 1/∆t divergences between different diagrams, such as for the sum of (2.28) and (2.39) in the present case. The problem is that there can be contributions coming from the pole at ∆t=0 that need to be accounted for. The previous paper [7] attempted to isolate these pole contributions by appropriate replacements of the form for each ∆t integral. The same method applied to the sequential diagrams would give a pole contribution to Φ + Ψ ≡ 1 2 [2 Re(xxyȳ + xxȳy) − IMC] + 2 Re(xyxȳ 2 ) of We will refer to this as a "1/π" contribution because of the single factor of π in the denominator. It turns out that both this result for the pole contribution to the sequential diagrams and the previous result [7] for the pole contribution to the crossed diagrams are incomplete: They miss some additional terms that involve 1/π 2 . The proper calculation of the pole terms is a lengthy enough issue that it would distract from our focus in this paper, and so here we will just quote results. The full discussion of why the previous analysis was incomplete, and of how to compute the full pole pieces (using dimensional regularization), is left to ref. [29]. With regard to sequential diagrams, the result is that (2.41) should be supplemented by the additional contribution For crossed diagrams, which are also needed for our total results of figs. 16 and 17, the corrected pole contribution will be given in section III B below. This correction is, in fact, critical to the Gunion-Bertsch-like cancellation of logarithms (1.10) discussed in Appendix B. 20

A. Sequential Diagrams
We now give a summary of our final formulas for the sequential diagrams of fig. 13, in the same style as the summary of the crossed-diagram contribution in section VIII of the preceding paper [7]. The two should be added together (along with the contributions involving 4-gluon vertices, like fig. 14, which we have left for the future).
The main result of this paper is ∆ dΓ dx dy sequential = A seq (x, y) + A seq (1−x−y, y) + A seq (x, 1−x−y) where A seq (x, y) is the result for (i) 2 Re(xyxȳ 2 ), for the large-N c color routing of fig. 24, plus half of (ii) 2 Re(xxyȳ + xxȳy) minus the corresponding idealized Monte Carlo result. That is, A seq corresponds to Ψ + Φ in the notation of (2.31) and (2.35). We will write this as D seq , defined as the integrand for xyxȳ 2 , is given by (ᾱ,β,γ), the (X seq , Y seq , Z seq ) and the I seq are given in appendix E by (E4), (E12), and (E13) respectively. The 4-particle evolution frequencies Ω ± are given by eq. (5.21) of ref. [7]. The M and other Ω here are and , Finally, the pole piece is In the last formulas, P (x) refers to the spin-averaged g→gg DGLAP splitting function 22 (3.10)

B. Crossed Diagrams
The formulas for crossed diagrams are summarized in section VIII of ref. [7], except that the pole contribution there needs to be corrected [29] by replacing the term 2 Re dΓ dx dy  in the formula (8.2) for A(x, y) in ref. [7] by The π/π 2 terms above are the same as in ref. [7] but are corrected by the addition of the 1/π 2 terms also included above.

IV. CONCLUSION
We have still not computed contributions from 4-gluon vertices such as fig. 14, which we expect will be important when no final gluon is soft, i.e. when y ∼ x ∼ z. (Conversely, we believe that these contributions will not be important when at least one final gluon is soft.) To use our results for energy loss calculations, we will also need to consider virtual corrections to single gluon bremsstrahlung. Both of these calculations have been left for future work.
One of the interesting aspects of our results in fig. 17 is the sign of the correction ∆ dΓ/dx dy due to overlapping formation times. We finish below with a qualitative picture of why the correction should be negative for y x, z. Unfortunately, we do not have any simple, compelling, qualitative argument for the other, more important case: How can we simply understand why the correction should be positive for y ∼ x z (and possibly also for y ∼ x ∼ z, depending on the size and sign of 4-gluon vertex contributions)?
A picture of why ∆ dΓ/dx dy is negative for y x, z Consider y x < z. As reviewed earlier, QCD formation times for softer emissions are smaller than formation times for harder emissions; so the formation time for emitting y is small compared to the formation time for emitting x. Now consider the case where these formation times overlap, as shown in fig. 25a. We've chosen to look at a case where the y emission happens some time δ after the midpoint of the x emission. The analogous contribution in the idealized Monte Carlo calculation is shown in fig. 25b. In the idealized Monte Carlo calculation, which ignores formation times, we have treated the "time" of the x and y emissions to be the midpoints of the corresponding formation times in fig. 25a. (This is the natural choice. Remember that it is also the choice we used in section II A 2 and technically justify in appendix D.) In the idealized Monte Carlo calculation, the chance of y emission from the x daughter is treated as completely independent from the chance of y emission from the z daughter, and x and overlapping formation times. t fx and t fy indicate the scale of the formation times for x and y emission respectively. (b) The corresponding idealized Monte Carlo contributions. δ is the time separation between the x emission and the y emission, defined in (a) as the separation between the midpoint times τ x and τ y of (2.15). As usual, transverse separations are highly exaggerated in the drawing. [There is no significance to (i) y being drawn angled up rather than down in the last term or (ii) the slightly exaggerated transverse separation of x in the drawing of (b) compared to (a). Both were done just to make the diagram clearer and less crowded.] these two probabilities are added together. But the separation of the x and z daughters is actually so small on these time scales that the y emission cannot resolve them as distinct particles when y x. Specifically, x emission is associated with transverse momenta of order Q ⊥,x ∼ qt form,x . Using the formation time (1.3), the corresponding scale of transverse separation of x from z during that formation time is b x ∼ 1/Q ⊥,x ∼ (xEq) −1/4 . When y x, this is indeed very small compared to the similar transverse scale characterizing the y emission during its formation time: b y ∼ (yEq) −1/4 . So, instead of seeing two daughter gluons (x and z), the y emission effectively sees a single adjoint-charge particle. (This resolution issue is similar to work by Mehtar-Tani, Salgado, and Konrad [30]. 23 ) On these time scales, there is therefore effectively only one particle providing a chance for initiating y emission rather than the two particles in the idealized Monte Carlo calculation associated with fig. 25b. That is, idealized Monte Carlo overcounts the probability for the y emission, and so the correction to idealized Monte Carlo should be negative for y x.   The following approximation reproduces the results of fig. 17 with a maximum absolute error 24 of 0.017 for all y > 10 −4 (assuming one permutes the final state gluons to choose y < x < z, just as in fig. 17): where the parameters each vary independently from 0 to 1. The numerical coefficients a mn and b mn are given in tables I and II. We have made no effort to make the approximation work well for y < 10 −4 .
2 Re(xyȳx + xȳyx) TABLE III: Coefficients of logarithms. The first column lists the subset of diagrams that are responsible for generating logarithms, using the notation of figs. 10-13 (to include permutations of x, y, and z). The second column identifies the UV-safe combinations that contain those diagrams, with A seq given by (3.2) of this paper. A is given by eq. (8.2) of ref. [7], as corrected by section III B of this paper. The last column is the corresponding leading-log contribution to ∆ dΓ/dx dy in units of

Appendix B: Logarithms and their cancellation
In this appendix, we discuss the behavior of individual contributions to ∆ dΓ/dx dy in the y x 1 limit. The individual contributions are logarithmically enhanced compared to the total, as in (1.10). In more detail, where ≈ indicates that we are only showing the leading logarithm. We can break the contributions down somewhat further, as shown in Table III. (One way to extract the coefficients in this table is from numerical extrapolation, as shown in fig. 26.)

DGLAP
There is a relatively simple way to understand some of these logarithms. Let's consider the 2 Re(xxyȳ+xxȳy) entry of Table III, which corresponds to the second and third diagrams of figs. 12 and 13. Since y is the softest high-energy particle in the process, it is the one most affected by the medium (which is, recall, why its formation time is the shortest). So focus attention on the contributions to this process where the only interaction with the medium occurs during the y emission. That is, consider fig. 27. In this picture, the x emission looks like an initial-state radiation correction to the underlying process of the y emission. Initialstate radiation generates a factorizable logarithm associated with DGLAP evolution. Let's calculate it. For the sake of familiarity, we'll start from the usual leading-order DGLAP evolution equation (though that is a somewhat backwards way to start given how DGLAP evolution is derived in the first place).  q E 1/2 ln x y . The results are plotted against 1/ log 10 (x/y), which allows for a simple linear extrapolation of the y → 0 limit for fixed x. Data points were calculated for x = 10 −3 , and we checked that results for x = 10 −4 would be indistinguishable to the eye. Note that the × and • data points are hard to distinguish because they are almost on top of each other. We have not shown points for A seq (y, x) + A seq (y, z) and A seq (z, y) + A seq (z, x): these would both lie right on top of the A seq (x, y) + A seq (x, z) points.
Consider the (vacuum) DGLAP evolution of an initial parton distribution function f (ζ, Q 2 ) for a parton of energy ζE. (We use the symbol ζ for momentum fraction here instead of the more traditional x or z just to avoid confusion with the specific use of x, y, and z ≡ 1−x−y in this paper as the momentum fractions of the three final-state gluons.) The basic DGLAP evolution equation takes the form where P is the relevant DGLAP splitting function. Now formally solve perturbatively by iteration, writing f = f 0 + f 1 + · · · with f 0 (ζ) = δ(1 − ζ) representing the case of no initial radiation. Then f 1 (corresponding to a single DGLAP emission from the initial state) is given by The solution is where we'll discuss what the reference momentum Q 2 0 should be in a moment. The result we want combines this probabilistic description of the initial x emission with the rate for the sub-process associated with the blob in fig. 27. That sub-process represents the rate of y emission in the medium, and so the leading-log estimate of the contribution to dΓ/dx dy associated with fig. 27 should be In the main text and in the preceding paper [7], our convention has usually been to identify the argument of the DGLAP splitting function with the momentum fraction x of the emitted gluon instead of with ζ = 1−x in fig. 27. Since the relevant splitting g → gg for this paper is symmetric, it doesn't matter in any case. So we'll henceforth rewrite P (1−x) above as P (x). The remaining job is to figure out the scales of the relevant upper and lower kinematic bounds on DGLAP virtuality Q 2 as applies to our problem, in order to know what ratio appears in the logarithm in (B5). In the high-energy, nearly-collinear limit, virtuality is related to off-shellness ∆E of energy by Q 2 ∼ E Q ∆E. In the case of fig. 27, E Q E for small x, where E is the initial parton energy. ∆E is related to time by the uncertainty principle, and so Q 2 ∼ E/∆t x , where ∆t x is the time scale associated with the x emission. So, (B5) can be translated to The time scale ∆t y of the underlying y emission is the formation time t form,y , which provides the lower limit on the time scale ∆t x for generating a DGLAP logarithm.
In the above DGLAP analysis based on fig. 27, we assumed that medium effects on the x emission could be ignored. This assumption is only valid for time scale small compared to the formation time t form,x for x emission. Another way to explain this infrared cut-off is to note that the logarithm in (B5) is a collinear logarithm. In vacuum, the x and ζ particle trajectories can be arbitrarily collinear, and there is no small-angle cut-off for collinear logarithms if we ignore the masses of particles. In medium, however, particles constantly change their direction over time, and so the x and ζ trajectories cannot remain collinear over arbitrarily large times. Formation times tell you how long you can wait before such changes destroy quantum coherence. The upshot is that the relevant range of ∆t for which a vacuum DGLAP analysis applies is t form,y ∆t x t form,x . (B7) Correspondingly the logarithm in (B5) is ln(t form,x /t form,y ). Combining t form,x ∝ √ x from (1.3) with the soft limits P (x) 2C A /x from (3.10) and dΓ dy from (2.24) and (2.25) then gives As promised, this agrees with the 2 Re(xxyȳ + xxȳy) entry of Table III. Readers may wonder why, in our discussion above (and especially in our discussion of Gunion-Bertsch cancellation below), we have focused on (i) vacuum x radiation from an underlying medium-induced soft y-emission process, rather than (ii) soft, vacuum y radiation from an underlying medium-induced x-emission process. The reason is that the latter contribution to real double bremsstrahlung dΓ/dx dy is sub-leading compared to the former by a factor of y/x. To understand this, consider the x↔y analog of (B6), which would contribute to dΓ/dx dy something of order α s P (y) dΓ dx × log. Since dΓ/dx ∼ α s /x 3/2 in units of q/E, this gives a contribution to dΓ/dx dy of order (α 2 s /x 3/2 y)×log, which is sub-leading compared to (B1) and (1.9).

Gunion-Bertsch
If we carry further the model we used in drawing fig. 27, we might expect the logarithmic corrections to the total rate to be similarly related to fig. 28. If we now abstract the y emission process and its attendant medium interaction as a net injection of momentum, then Fig. 28 looks perhaps analogous to the process of fig. 29. The latter figure shows the type of non-Abelian radiation process long ago considered by Gunion and Bertsch [25] in the context of 2-particle collisions in vacuum. In the high-energy, nearly-collinear limit, the amplitude takes the form 25 (B10) where ⊥ is defined relative to the initial particle direction, l ⊥ and a characterize the transverse momentum transfer and adjoint color index associated with the collision, and k ⊥ and b characterize the transverse momentum and adjoint color index of the final bremsstrahlung 25 Our k ⊥ is the q ⊥ of Ref. [25]. We have not bothered to show the matrix element representing the cross in our fig. 29, which can be factorized out. gluon. In the k ⊥ l ⊥ limit, if we keep only the leading behavior of each individual diagram, the corresponding rate is proportional to (B11) Any individual term, such as generates a factor proportional to g 2 ln(k max ⊥ /k min ⊥ ), which is the same sort of logarithmic factor found in (B5). However, all the terms in the amplitude in (B11) cancel each other, and so all of these logarithmic factors coming from k ⊥ l ⊥ cancel each other in the total rate. Moreover, if the initial particle is a gluon, so that the generators T above are adjointrepresentation generators, then the squares of each of the individual three terms in (B11) are respectively proportional to which (for the adjoint representation) are all the same! In our analogy, this corresponds to the equality of the last three rows of Table III. Similarly, the cross-terms are proportional to which all equal the negative of (B13), analogous to the first three rows of Table III. So the cancellation of the logarithms in Table III appears to be of the Gunion-Bertsch type. What about the condition k ⊥ l ⊥ for Gunion-Bertsch cancellation? First consider that for k ⊥ l ⊥ and small x, the off-shellness ∆E of energy associated with the intermediate lines in fig. 28 are all of order k 2 ⊥ /xE, which corresponds to a time scale ∆t x ∼ xE/k 2 ⊥ . The condition k ⊥ l ⊥ is therefore equivalent to the condition ∆t x xE/l 2 ⊥ . In the context of fig. 28, l ⊥ consists of the combination of the transverse momentum carried away by the y bremsstrahlung and the transverse momentum ∼ qt form,y injected by the medium during the y formation, which are parametrically the same and give l ⊥ ∼ qt form,y . So the condition k ⊥ l ⊥ is equivalent to Using the formation time (1.3), this may be recast as Given our original assumption that y x, the above condition is automatically satisfied for the range (B7) of values that generated the logarithm in our case. The moral of this story is that the logarithms characterized by Table III arise from a kinematic regime that is analogous to the cancellation of k ⊥ l ⊥ logarithms in Gunion-Bertsch.

Independent emission model
There is an alternative picture for understanding and interpreting two of the entries in Table III, which we offer here as a complement to the previous discussion.

a. The approximation
Imagine that we tried to estimate double bremsstrahlung in the y x 1 limit by assuming that the x and y bremsstrahlung amplitudes were completely independent from each other and both given simply by single-bremsstrahlung formulas. This approximation would be somewhat similar to the small x and y limit of the idealized Monte Carlo approximation (2.8) except that we don't care about the ordering of τ x and τ y : We will find that this is a useful approximation for the processes indicated in the first two rows of Table III, for which the time integrations above must be correspondingly restricted.  between the x emissions. This time ordering modifies (B17) to Why have we written •yȳ• + •ȳy• in the subscript above instead of, say, xyȳx + xȳyx? In our limit y x 1, we have z ≡ 1−x−y 1−x. For single bremsstrahlung, there is no difference between referring to the corresponding diagram ( fig. 21) as xx or zz (with z ≡ 1−x in the single bremsstrahlung context). We could use either notation to describe the x emission in the independent bremsstrahlung approximation we have outlined above. When we talk about actual double bremsstrahlung diagrams, however, there is a difference. Let's look at those diagrams: The interference processes listed in the first two rows of Table III can be drawn in the form of fig. 30. However, for the reasons given in section IV, we expect that the y emission cannot resolve the two daughters (x and 1−x) of the x emission process when y x 1. So y emission from those two daughters is approximately equivalent to y emission from a single particle of energy E (and so is similar to the cases where y emission occurs instead from the initial particle of energy E). That is, we can redraw the sum of diagrams in fig. 30 as fig. 31, and it is this sum that should correspond to the independent approximation given by the right-hand side of (B18). The subscript notation on the lefthand side of (B18) means that this approximation corresponds to the diagrams of fig. 30 with "•" summed over the cases • = x and • = z.  fig. 30 in the limit y x 1 where the y emission cannot resolve the x and 1−x daughters of the x emission process and couples to them (magenta ovals) as to a single adjoint-representation particle.

b. Evaluation and Comparison
Comparison of (2.6) and (2.24) identifies Using (3.10), the independent emission contribution (B18) for y x 1 is then Re Ω 2 y csc 2 (Ω y ∆t y ) Re ln where we define the small x, y expressions The integrand falls exponentially for ∆t y Ω −1 y , and so Ω x ∆t y 1 in the y x 1 limit, giving This integral has a ∆t→0 divergence. In order to regulate divergences of individual diagrams, we have routinely subtracted out the vacuum contribution to those diagrams. Here, we will subtract out the vacuum contribution to the independent y emission, replacing (B22) by dΓ dx dy independent 2 Re(•yȳ•+•ȳy•) where γ E is the Euler-Mascheroni constant. (We leave a more careful and convincing treatment of the ∆t→0 regularization to other work [29].) Using (B21), the above result can be recast as The coefficient −1 of the logarithm indeed matches up with the sum of the first two rows in Table III, as promised. We should mention that it is also possible to extract the integral expression in (B23) not just from the simple independent bremsstrahlung approximation above but also from the full integral expressions for the general result for 2 Re(xyȳx + xȳyx) and 2 Re(zyȳz + zȳyz) given in the preceding work [7]. This can be done by analytically extracting the behavior of those expressions in the limit y x 1 while formally treating the integration variable ∆t y ∼ Ω −1 y ∝ √ y. We will not give details here. Individually, 2 Re(xyȳx + xȳyx) and 2 Re(zyȳz + zȳyz) each contribute half of the limiting result (B24) in this analysis.
In the first row of Table III, we did not explicitly list the 2 Re(xȳxy) contribution to the group of diagrams represented by A(x, y). 26 To understand this omission, note that the origin of the logarithm in (B23) can be traced back to the ∆t x − ∆t y ∆t x factor inside the last integral in (B18). This factor arises because, for •yȳ• + •ȳy•, the short-duration y emission can happen at any time in the middle of the longer-duration x emission. For the xȳxy process shown in fig. 10, however, (and similarly for zȳzy) the short-duration y emission is forced by the time-ordering of the diagram to occur only at the very end of the x emission interval, and so there is no comparable factor and so no logarithmic enhancement.

c. Other processes
One may investigate similar approximations of other processes in Table III. However, we do not know of any sensible way to apply the independent bremsstrahlung approximation to the third row of the table, which refers to 2 Re(xzzx + xzxz + zxxz + zxzx). These diagrams are depicted in fig. 32 in a fashion similar to fig. 30. These are interference diagrams that combine with the Monte Carlo related diagrams of fig. 33 to partly suppress the result along the lines discussed in section IV. Roughly speaking, the y emission's inability to individually resolve the x and 1−x daughters in the y x limit partly suppresses the rate. This combined effect is depicted in fig. 34. The problem with applying an independent bremsstrahlung approximation to the picture in fig. 34, however, is that there is nothing in that approximation that would constrain how large could be the separation of the y emission from the x emission. Consider the two daughters of the x emission. The issue of how y emission from one daughter decoheres with y emission from the other, for times t form,x after the last x emission, is a complicated matter that is simply not captured by the independent bremsstrahlung approximation that was used so successfully above for fig In this appendix, we show how to start with the methods and diagrammatic rules of the preceding paper [7] and obtain the result for the xxyȳ interference given by (2.18).
Using the notation of ref. [7], the xxyȳ interference of fig. 35 gives |i δH|Bȳ Bȳ, tȳ|B y , t y B y |−i δH| for the differential probability. Here, all transverse positions (B x , Bx, B y , Bȳ) are implicitly integrated over, and δH is the part of the full Hamiltonian of the theory that contains the splitting vertices for high-energy particles. One non-trivial normalization factor above has to do with the intermediate 2-particle (effectively 0-particle) state | for tx < t < t y . This state was normalized in the preceding paper [7] as where, in this context, (x 1 , x 2 ) = (x 1 , −x 1 ) are the longitudinal momentum fractions of the two particles in the 2-particle state. This normalization was both natural and convenient in the discussion of ref. [7], and it is similarly convenient here, for the normalization of the initial and final states in (C1). However, having a non-trivial normalization means that the correct way to project out the 2-particle (effectively 0-particle) intermediate state, as in (C1). The effectively 0-particle intermediate state | does not have any time dependence in our formalism, and so the integrand in (C1) does not care how far apart in time tx and t y are. That is, we have not bothered to write a factor of 1 in (C1) representing the time evolution of this intermediate state: This time independence is consistent with interpreting this diagram as representing two consecutive splittings that are completely independent from each other. Using the rules developed in the preceding paper [7], fig. 35 and (C1) become where P is defined in terms of DGLAP splitting functions in section IV.E of ref. [7]. The overall (1 − x) −2 factor above is the intermediate-state normalization factor | −1 from (C1).
In the application in this paper, the color representation R of the solid line in fig. 35 will be adjoint, and so C R above is C A . Now turn to the helicity sums, which we will relate to spin-averaged DGLAP splitting functions. To this end, it is useful to first use rotational invariance in the transverse plane to write and so Up to normalization factors in the definition of P, the two factors above each represent a spin-averaged DGLAP splitting function. Taking the precise definition of P from the preceding paper [7], we specifically obtain Substituting (C11) for (C8) in (C7) then yields 27 with y defined by (2.1). Now using the definition (2.7) of dΓ/dx d∆t, (C12) will give (once we clarify below some issues about normalization) which is the result claimed in (2.18) in the main text (after using time translation invariance to eliminate one of the time integrals and so change probability distribution dI/dx dy to rate distribution dΓ/dx dy). Let's now address the normalization issue, since (C12) and (C13) may not look like they exactly match up. In particular, one might think that the formula for [dΓ/dy ∆t] (1−x)E is obtained from the definition (2.7) of [dΓ/dx d∆t] E by simply substituting (1 − x)E, y for (E, x). But this overlooks the fact that the effective 1-particle quantum mechanics variables B here and in the previous paper [7] are defined using the longitudinal momentum fractions defined with respect to the initial parent E of the full double-splitting process. In the case of the t y < t < tȳ evolution in fig. 35 and (C5), 27 An alternative way one could get to this result from (C7) is to use (E3) to write (C8) as 4(ᾱ + 1 2β + 1 2γ ). Then use (E5) to arrive at (C12).
where b z and b y are the transverse positions of the z and y daughters. In contrast, the B in the corresponding single-splitting calculation would be which is (other than the choice of notation) the same thing as b z −b y . The correct translation of (2.7) is then (C16) Recalling that the states |B are normalized so that B|B = δ This indeed matches up with the corresponding factor in (C12) to give (C13).
Appendix D: Justification of the prescription (2.15) We now outline how to justify the midpoint prescription (2.15). As discussed in the main text, the starting point is to imagine a medium with, for example, In the limit of large T,q is approximately constant over any formation time, and so we may still analyze the effects of overlapping formation times in the constant-q approximation used for detailed calculations in this paper. But, at the same time, (D1) helps by providing a consistent large-time regulator for independent splittings that are far separated in time.
There is a subtlety, however, to approximatingq as constant over formation times. Consider, formally, the contribution of xxyȳ to the total probability I for double splitting: When evaluating dΓ/dy d∆t y in the integrand, should we takeq to beq(t y ) orq(tȳ) or pick some other time in between? [A similar question arises for dΓ/dx d∆t x .] This may at first sound unimportant because the integrand falls off exponentially when ∆t y exceeds the formation time, and so (except for negligible corrections) t y and tȳ can be at most roughly a formation time apart. That meansq(t y ) →q(tȳ) as T → ∞. The difference between q(t y ) andq(tȳ), and so the different values we might get for dΓ/dy d∆t y , scale at worst like T −1 . The problem is that this correction can accumulate additively when we integrate over time. For simplicity of argument, let's focus on the integral over t y while holding t x , tx, and ∆t y ≡ tȳ − t y fixed. If we choose to evaluateq at t y , then the dt y integral has the functional form for some function f (q). If we instead takeq(tȳ), it is The difference between these two prescriptions is which does not vanish as T → ∞, and so the difference of the prescriptions forq cannot be ignored. So which choice is correct? We could avoid this question by computing dΓ/dy d∆t y [and similarly dΓ/dx d∆t x ] without making the approximation thatq is constant over ∆t, but that makes the calculation unnecessarily difficult. Instead, consider the following. For simplicity of argument, focus on a region of time whereq(t) is monotonically decreasing, i.e. t > 0 in (D1). The problem with the choiceq(t y ) for the time interval (t y , tȳ) is then that it overestimatesq(t) for the entire interval. The problem with the choiceq(tȳ) is that it underestimatesq(t) for the entire interval. However, if we choose the mid-point,q(τ y ) with τ y ≡ (t y + tȳ)/2, we will in almost equal measure overestimateq in the first half of the interval and underestimate it in the second half. That is, the error we make in dΓ/dy d∆t y will be O(T −2 ) instead of O(T −1 ). 28 This is small enough that the accumulated error after integrating t y over time will not affect our final results in the T → ∞ limit.
The upshot is that, in the xxyȳ calculation,q should be evaluated at the midpoint times τ x and τ y . In the single-splitting calculation, it should also be evaluated at the midpoint time, and so the formula used for the idealized Monte Carlo calculation also corresponds to usingq(τ x ) andq(τ y ), where τ x and τ y are the emission times in the idealized Monte Carlo. We could now go through all the steps from (2.17) to (2.23) 29 explicitly specifying the use 28 For a more explicit discussion of evaluating single-splitting rates in the case of generic, time-varyingq, see ref. [33]. In terms of explicit formulas, the argument about error sizes used above can be reproduced from eq. (3.8) of ref. [33] by expanding ω 2 0 (t) = ω 2 0 (t mid ) + (t − t mid ), where ∼ T −1 is small and t mid ≡ (t 1 + t 2 )/2. Note that the integral and integrand of eq. (3.8) of ref. [33] are symmetric under t 1 ↔ t 2 by eq. (3.24) of that reference. By the corresponding anti-symmetry of the (t − t mid ) term in the expansion of ω 2 0 (t) under reflection about t mid , there can be no O( ) contribution to the probability. The first correction will instead be O( 2 ). 29 We are glossing over a few details here, having conflated calculations of contributions to the doublesplitting probability (D2) with calculations of contributions to the double-splitting rate and thence ∆ dΓ/dx dy (2.23). Giving the time dependence (D1), one should initially discuss computing probabilities rather than rates. However, once one computes the correction ∆ dI/dx dy to the probability due to overlapping formation times, one may then define the corresponding correction ∆ dΓ/dx dy by casting the probability into the form ∆ dI/dx dy = dt ∆ dΓ/dx dy in the case of generic, arbitrarily slowly varying choices ofq(t) (that fall off to zero as t → ±∞). When thinking about the size of effects that must be FIG. 36: Our labeling conventions for the xyxȳ 2 interference diagram. (The only difference between xyxȳ 2 and xyxȳ 1 is the order (x 1 , x 2 , x 3 , x 4 ) we have chosen for labeling the particles during the 4-particle evolution t y < t < tx.) ofq(τ x ) andq(τ y ). However, once we get to the final result (2.23), which is finite as T → ∞, we no longer have to worry about this distinction. We can then simply take T → ∞ and callq a constant.
The color factor 1 2 C 2 A in front is half of the naive color factor d −1 fig. 36 because xyxȳ 2 given by (E1) represents the contribution from only one of the two large-N c color routings of xyxȳ.
The helicity sums are not as simple as in the xxyȳ case of appendix C because we do not have any useful analog of (C6) to start the simplification. However, we may again average over the initial helicity h i , and rotational symmetry in the transverse plane implies that the initial helicity average of hx,hy,hz h above must have the form α(x, y) δn n δm m +β(x, y) δnmδ nm +γ(x, y) δn m δ nm , for some functionsᾱ,β, andγ, analogous to eq. (4.38) of ref. [7]. Using the explicit formulas for the P's for the case where all high-energy particles are gluons, 30 It will be useful to note for future reference that a certain combination of (ᾱ,β,γ) can be written in terms of spin-averaged DGLAP splitting functions: .
However, they are symmetric under y ↔ z ≡ 1−x−y.
In the harmonic oscillator approximation, we may do the first and last time integrals (t x and tȳ) as in section V.A of ref. [7], with result (ᾱδn n δm m +βδnmδ nm +γδn m δ nm ) where ∆t ≡ tx − t y and Ω i and Ω seq f are respectively the Ω E,x and Ω (1−x)E,y of (2.25) and (2.27). We put the superscript "seq" on M seq f and Ω seq f to indicate that they are different than the corresponding M f and Ω f in the calculation of the canonical diagram xyȳx in ref. [7].
Unlike in the xyȳx analysis of the preceding paper [7], the relevant variables (C 41 , C 23 ) are the same on both sides of Cx 41 , Cx 23 , tx|C y 41 , C y 23 , t y above. In order to keep the notation similar to that of the preceding paper, it is useful to consider the variables ordered on the two ends as Cx 23 Cx 41 and C y 41 C y 23 (E9) so that the derivatives in (E7) hit the lower element of each pair (as in that paper). The appropriate transformation matrix between these variables and the normal modes for the 4-particle evolution is at t y , in the notation of the preceding paper. The explicit formula is given by eqs. (5.24,5.28,5.33) of that paper. Given the choice (E9), the corresponding transformation at tx is simply The rest of the analysis proceeds as in section V.C of ref. [7] but using in result (2.36) of this paper. 32 The I seq 's in (2.36) are defined the same way as the I's of ref. [7] but using the appropriate (X seq , Y seq , Z seq ) here instead of (X, Y, Z) there: . (E13e) The 4-particle evolution frequencies Ω ± appearing in (2.36) are given by eq. (5.21) of ref. [7].
In the current context, sgn M i = +1 and sgn M seq f = +1. We will not need to make use of this result in any other context, and so we will drop these signs. Using the above expansions in (2.36) and (E13) then gives which may be recast into the form (2.37) with the aid of (2.33), (E5) and (E8).
Appendix F: Transverse momentum integration vs. Ref. [22] Consider the case of independent single-splitting processes-what we've called the idealized Monte Carlo (IMC) picture. For democratic branchings (no daughter soft), there is a simple qualitative argument for the IMC picture: The chance of bremsstrahlung is roughly α s per formation time t form , and so the typical time between (non-soft) splittings is τ rad ∼ t form /α s , which is parametrically large compared to the formation time t form when α s at the relevant energy scale is small. For the LPM effect, the formation time represents the time over which the splitting process is quantum mechanically coherent, and so events which are separated by more than that time should be quantum mechanically independent. The effects we calculate in this paper, in contrast, are suppressed by a factor of t form /τ rad ∼ α s , which is the chance for two consecutive splittings to overlap.
In the context of QCD, the basic scales described above have been known since the pioneering work of Baier et al. [31] and Zakharov [32] (BDMPS-Z). There have been many works on extending (in various limits) the BDMPS-Z analysis to include transverse momentum distributions, but the issue we want to mention here arises in the work of Blaizot et al. [22]. One of the issues that concerned them was how long color coherence between the daughters might survive after a splitting, and how formally to show that it does not interfere with the quantum-mechanical independence of splittings. They found that color decoherence occurs over a time of order t form after the time that we call tx in our fig. 21, and so all is well qualitatively concerning the separation of scales justifying the leadingorder approximation that splittings are independent. The analysis of this decoherence was somewhat complicated, even in large N c .
However, as we'll now discuss, this late-time color decoherence issue (and so our need to calculate it) disappears completely if one integrates the splitting rate over the transverse momenta of the daughters after the splitting, which means that we can avoid calculating the details of late-time color decoherence. Let's start by discussing the issue in the language used by ref. [7], which developed many of the methods used in our paper here. The issue is treated in section IV.A of ref. [7] in the context of double splitting, where it is shown using a unitarity argument that, provided one integrates over the final transverse momenta of the daughters, then one may ignore what happens to a daughter after it has been emitted in both the amplitude and the conjugate amplitude. So, for instance, consider the xyȳx diagram of figs. 10 and 11. For this diagram, one may ignore what happens to the y daughter after the corresponding splitting time (tȳ) in the conjugate amplitude, and one may additionally ignore what happens to the other two daughters after the corresponding splitting time (tx) in the conjugate amplitude. The same arguments apply to single splitting, and one may ignore what happens to the two daughters after the time tx in the xx diagram of fig. 21. If we did not integrate over transverse momenta, then we would have to worry about later interactions, such as shown in fig. 37 for the case of single splitting, which represents the same interactions as fig. 9 of Blaizot et al. [22]. The development of the system for times later than tx in our figure is what Blaizot et al. refer to as "Region III" in their discussion, and the tricky color coherence issues are associated with the time interval tx < t < t 3 , which they refer to as the "non-factorizable" piece of the calculation.
In the formalism of Blaizot et al., what happens when we do integrate over final transverse momenta? Integrating over final transverse momenta k a and k b in their differential crosssection (3.13) is equivalent, in the transverse position space used in their (4.1), to equating (ȳ a ,ȳ b ) = (y a , y b ) and then integrating over (y a , y b ). In their formula (4.15) for the effects of late-time correlations such as shown in our fig. 37, the factor (y a ;ȳ a |S (2) (t L , t 3 )|z a ;z a ) then forces z a =z a . (Similarly, z b =z b .) This makes the T (Z) of their (4.16) vanish, which means that their (4.15) vanishes: what they call the "non-factorizable" contributions vanish when integrated over final transverse momenta.