The LPM effect in sequential bremsstrahlung: nearly complete results for QCD

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. We continue study 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-emission approximations. Previous work has computed overlap effects for double splitting $g \to gg \to ggg$. To make use of those results, one also needs calculations of related virtual loop corrections to single splitting $g \to gg$ in order to cancel severe (power-law) infrared (IR) divergences. This paper provides calculations of nearly all such processes involving gluons and discusses how to organize the results to demonstrate the cancellation. In the soft emission limit, our results reproduce the known double-log behavior of earlier authors who worked in leading-log approximation. We also present a first (albeit numerical and not yet analytic) investigation of sub-leading, single IR logarithms. Ultraviolet divergences appearing in our calculations correctly renormalize the coupling $\alpha_{\rm s}$ in the usual LPM result for leading-order $g \to gg$.

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. The goal of this paper is to (i) present nearly complete results for the case of two overlapping gluon splittings (e.g. g → gg → ggg) and (ii) confirm that earlier leading-log results for these effects [7][8][9] are reproduced by our more-complete results in the appropriate soft limit. As a necessary step, we discuss how to combine the effects of overlapping real double splitting (g → gg → ggg) with corresponding virtual corrections to single splitting (e.g. g → gg * → ggg * → gg) to cancel spurious infrared (IR) divergences. In our analysis of virtual corrections, we will also verify that we reproduce the correct ultraviolet (UV) renormalization and running of the QCD coupling α s associated with the high-energy vertex for single splitting.
In this paper, we will present the formulas for the building blocks just discussed, but we leave application of those formulas to later work. In particular, one of the ultimate motivations [10] of our study is to eventually investigate whether the size of overlap effects is small enough to justify a picture of parton showers, inside a quark-gluon plasma, as composed of individual high-energy partons; or whether the splitting of high-energy partons is so strongly-coupled that high-energy partons lose their individual identity, similar to gaugegravity duality studies [11][12][13][14] of energy loss. But, as will be discussed in our conclusion, further work will be needed to answer that question.
As a technical matter, our calculations are organized [15] using Light-Cone Perturbation Theory (LCPT) [16][17][18]. 2 As we will explain below, the "nearly" in our claim of "nearly complete results" refers to the fact that we have not yet calculated, for QCD, contributions from diagrams that involve "instantaneous" interactions in Light-Cone Perturbation Theory. The effects of such diagrams have been numerically small in earlier studies of overlap effects in QED [15], and they do not contribute to our check that our results agree with earlier leading-log calculations. For these reasons, and because analysis of the non-instantaneous diagrams is already complicated, we leave the calculation of instantaneous diagrams for QCD to later work. For similar reasons, we also leave to later work the effect of diagrams involving 4-gluon vertices, like those computed for real double gluon splitting in ref. [20].
We make a number of simplifying assumptions also made in the sequence of earlier papers [15,[21][22][23] leading up to this work: We take the large-N c limit, assume that the medium is thick compared to formation lengths, and use the multiple-scattering (q) approximation appropriate to elastic scattering of high-energy partons from the (thick) medium. All of these simplifications could be relaxed in the context of the underlying formalism used for 1 The papers of Landau and Pomeranchuk [1] are also available in English translation [3]. The generalization to QCD was originally carried out by Baier, Dokshitzer, Mueller, Peigne, and Schiff [4,5] and by Zakharov [6] (BDMPS-Z). 2 For readers not familiar with time-ordered LCPT who would like the simplest possible example of how it reassuringly reproduces the results of ordinary Feynman diagram calculations, we recommend section 1.4.1 of Kovchegov and Levin's monograph [19]. calculations, 3 but practical calculations would then be quite considerably harder; so we focus on the simplest situation here.
A. The diagrams we compute Previous work [21][22][23] has computed overlap effects for real double gluon splitting (g → gg → ggg) depicted by the interference diagrams of figs. 1 and 2. Each diagram is timeordered from left to right and has the following interpretation: The blue (upper) part of the diagram represents a contribution to the amplitude for g → ggg, the red (lower) part represents a contribution to the conjugate amplitude, and the two together represent a particular contribution to the rate. Only high-energy particle lines are shown explicitly, but each such line is implicitly summed over an arbitrary number of interactions with the medium, and the diagram is averaged over the statistical fluctuations of the medium. See ref. [21] for details. For real double gluon splitting, we will refer to the longitudinal momentum fractions of the three final-state gluons as x, y, and z ≡ 1−x−y (1.1) 3 In particular, for a discussion of how one could in principle eliminate the large-N c approximation, see refs. [24,25]. Time-ordered diagrams for the leading-order rate for single gluon splitting. [22].
relative to the initial gluon. Also, our nomenclature is that figs. 1 and 2 are respectively called "crossed" and "sequential" diagrams because of the way they are drawn. For the case of sequential diagrams ( fig. 2), it is possible for the two consecutive splittings to be arbitrarily far separated in time, in which case their formation times do not overlap. The effect of overlapping formation times in this case is then determined by subtracting from the sequential diagrams the corresponding results one would have gotten by treating the two splittings as independent splittings. Details are given in ref. [22], along with discussion of physical interpretation and application. 4 Whenever such a subtraction needs to be made on a double-splitting differential rate dΓ, we will use the symbol ∆ dΓ to refer to the subtracted version that isolates the effect of overlapping formation times.
In the limit that one of the three final-state gluons-say y-is soft, it was found [22] that the overlap effect on real double splitting behaves parametrically as 5 As we'll review later, the y −3/2 behavior would lead to power-law infrared divergences in energy loss calculations. Very crudely analogous to what happens in vacuum bremsstrahlung in QED, where there are (logarithmic) infrared divergences that cancel in inclusive calculations between real and virtual emissions, we need to supplement the real double emission processes (g → ggg) by a calculation of corresponding virtual corrections to the single emission process (g → gg) of fig. 3. The virtual processes that we calculate in this paper are shown in fig. 4 (which we call Class I) and fig. 5 (which we call Class II). There are also cousins of the Class I diagrams generated by swapping the two final state gluons (x → 1−x), two examples of which are shown in fig. 6. For Class II diagrams, such a swap does not generate a new diagram. In total, these sets of virtual diagrams include all one-loop virtual corrections to single splitting except for processes involving instantaneous interactions or fundamental 4-gluon vertices. As mentioned previously, we leave the latter for future work. A few examples are shown in fig. 7. The "instantaneous" interactions (indicated by a propagator crossed by a bar) are instantaneous in light-cone time and correspond to the exchange of a longitudinallypolarized gluon in light-cone gauge. See ref. [15] for examples of such diagrams evaluated in QED.
We should clarify that, physically, the power-law divergences of (1.2) as y→0 are not actually infinite. The scaling (1.2) depends on theq approximation, which breaks down when the soft gluon energy yE becomes as small as the plasma temperature T . 6 In the high-energy limit, however, the cancellation of such power-law contributions to shower development, even if only a cancellation of contributions that are parametrically large in energy rather than truly infinite, will be critical to extracting the relevant physics that survives after the cancellation. In this paper, we will be able to ignore the far-infrared physics (meaning scale T ≪ E) that regulates the power-law divergences and can simply analyze the cancellation of power-law divergences in the context of theq approximation appropriate for the high-energy behavior. Longitudinal gluon exchange is represented by a vertical (i.e. instantaneous) line that is crossed by a black bar, following the diagrammatic notation of Light-Cone Perturbation Theory.

B. Infrared Divergences
We will later discuss the calculation of the differential rates associated respectively with the real double emission diagrams of fig. 1 plus fig. 2, the Class I virtual correction diagrams of fig. 4, and the Class II virtual correction diagrams of fig. 5. But here we first preview some results concerning infrared divergences.
In the virtual diagrams of figs. 4 and 5, the virtual loop longitudinal momentum fraction y in the amplitude or conjugate amplitude needs to be integrated over, and it will be convenient to introduce the notation [dΓ/dx dy] virt I and [dΓ/dx dy] virt II for the corresponding integrands of that y integration. Our calculations are performed in Light Cone Perturbation Theory, in which every particle line (virtual as well as real) is restricted to positive longitudinal momentum fraction. The structure of the Class I diagrams of fig. 4 then forces 0 < y < 1−x, whereas the structure of the Class II diagrams of fig. 5 forces 0 < y < 1 instead. So, in our notation, We will later give detailed discussion of how infrared divergences appear in various calculations associated with shower development, but a good starting point is to consider the net rate [dΓ/dx] net at which all of the processes represented by figs. 1-6 produce one daughter of energy xE (plus any other daughters) from a particle of energy E, for a given x. That's given by where the first term is the rate of the leading-order (LO) g → gg process of fig. 3, and where the next-to-leading-order (NLO) contribution is 7 (1.5b) [See appendix B for more discussion.] The bars above LO and NLO in (1.5) are a technical distinction that will be discussed later and can be ignored for now. In the integrals above, some virtual or final particle has zero energy at both the lower and upper limits of the y integrations, and so both limits are associated with infrared divergences. In order to see how divergences behave, it is convenient to use symmetries and/or change of integration variables to rewrite the integrals so that the infrared divergences of [dΓ/dx] NLO net are associated only with y → 0 (for fixed non-zero x < 1). In particular, ( where contributions from virtual and real double splitting processes appear in the respective combinations The θ(· · · ) in (1.7) represent unit step functions [θ(true) = 1 and θ(false) = 0], and they just implement upper limits on the y integration. The advantage of using the θ functions is so that we can combine all the integrals: the integrals for the separate terms each have power-law IR divergences, but whether or not those divergences cancel is now just a question of the y → 0 behavior of the combined integrand of (1.7).
In the limit y → 0 for fixed x, the integrand of (1.7) approaches Using the symmetry of the g → ggg rate (1.8b) under permutations of x, y, and z = 1−x−y, we have r(x, y) = r(1−x−y, y) ≃ r(1−x, y) for small y, and so (1.9) approaches v(x, y) (1.10) By (1.2), r(x, y) ∼ y −3/2 for small y, and so the integral of r(x, y) in (1.7) has a power-law IR divergence proportional to 0 dy/y 3/2 . From the full results for rates that we calculate in this paper, we find that the y −3/2 behavior cancels in the combination v(x, y) + 1 2 r(x, y) appearing in (1.10). We also find that left behind after this cancellation is, at leading logarithmic order, v(x, y) which generates an IR double log divergence when integrated over y. As we discuss later, this result, applied to (1.7), exactly matches leading-log results derived earlier in the literature [7][8][9] and so provides a crucial check of our calculations. Though it should be possible to extract (1.11) from our results analytically, so far we have only checked numerically. 8 Fig. 8   vs. ln y for a sample value of x. According to (1.11), the slope of (1.12) vs. ln y should approach −1 as ln y → −∞, which we show in fig. 8 by comparison to the straight line. We hope in the future to also provide exact analytic results for single-log divergences that are subleading to the double-log divergence. For now we only have numerical results for those, which we present later with an examination of how well those numerical results fit an educated guess for their analytic form.

C. Outline
The new diagrams needed for this paper are the virtual diagrams of figs. 4 and 5. In the next section, we discuss how we can avoid calculating any of these diagrams from scratch. All of the g → gg QCD virtual diagrams can be obtained by either (i) transformation from known results for the g → ggg QCD diagrams of figs. 1 and 2 or (ii) by adapting the known result for one QED virtual diagram.
In section III, we go into much more detail about how to organize IR divergences in calculations related to energy loss. We also show that the double-log behavior (1.11) is equivalent to earlier leading-log results.
Section IV presents numerical results for sub-leading single-log divergences and shows that the numerics fit very well, but not quite perfectly, a form one might guess based on the physics of double-log divergences.
The formalism and calculations that have led to our results for rates have spanned many papers, and one can reasonably worry about the possibility of error somewhere along the way. Section V provides a compendium of several non-trivial cross-checks of our results. Section VI offers our conclusion and our outlook for what needs to be done in future work. Appendix A contains a complete summary of all our final formulas for rates. Many technical issues, derivations, and side investigations are left for the other appendices.

A. Symmetry Factor Conventions
Before discussing how to find formulas for differential rates, we should clarify some conventions. Note each virtual diagram in fig. 5, as well as the second row of fig. 4, has a loop in the amplitude (an all-blue loop) or conjugate amplitude (an all-red loop) that should be associated with a diagrammatic loop symmetry factor of 1 2 . Our convention in this paper is that any such diagrammatic symmetry factor associated with an internal loop is already included in the formula for what we call ∆ dΓ/dx dy in (1.4). Note that the loops in the first row of fig. 4 do not have an associated symmetry factor.
In contrast, we do not include any identical-particle final-state symmetry factors in our formulas for differential rates. These must be included by hand whenever integrating over the longitudinal momentum fractions of daughters if the integration region double-counts final states. For example, the total rate for real double-splitting g → ggg is formally given by because the integration region used above covers all 3! permutations of possible momentum fractions x, y, and z = 1−x−y for the three daughter gluons. Similarly, for g → gg processes, formally We use the caveat "formally" because the total splitting rates Γ and ∆Γ above are infrared divergent, but they provide simple examples for explaining our conventions.

B. Relating virtual diagrams to previous work
In the context of (large-N f ) QED, ref. [15] showed how many diagrams needed for virtual corrections to single splitting could be obtained from results for real double splitting via what were named back-end and front-end transformations. For the current context of QCD, figs. 9 and 10 depict diagrammatically how all but two of the Class I and II virtual diagrams we need (figs. 4 and 5) can be related to known results for crossed and sequential g → ggg diagrams (figs. 1 and 2) using back-end and front-end transformations, sometimes accompanied by switching the variable names x and y and/or complex conjugation. Diagrammatically, a back-end transformation corresponds to taking the latest-time splitting vertex in one of our rate diagrams and sliding it around the back end of the diagram from the amplitude to the conjugate-amplitude or vice versa. Diagrammatically, a front-end transformation The black arrows indicate moving the latest-time (or earliest-time) vertex using a back-end (or front-end) transformation [15].
corresponds to taking the earliest-time splitting vertex and sliding it around the front end of the diagram.
In terms of formulas, the only effect of a back-end transformation is to introduce an overall minus sign in the corresponding formula for dΓ/dx dy [15] requiring the longitudinal momentum fractions of the lines of the diagrams to match up requires replacing where E is the energy of the initial particle in the real or virtual double-splitting process. See section 4.2 of ref. [15] for a more detailed discussion. There is also an overall normalization factor associated with the transformation that, for our case here where all the particles are gluons, amounts to 9 dΓ dx dy in 4−ǫ spacetime dimensions. The overall factor (1−x) −ǫ will be relevant because we will use dimensional regularization to handle and renormalize UV divergences in our calculation. We should note that there are a few additional subtleties in practically implementing frontend transformations, which we leave to appendix D. As an example of (2.7), the relation depicted by the first case of fig. 10 gives See appendix H of ref. [15], especially eqs. (H.13) and (H.14) there. In (H.13) of ref. [15] there was additionally an overall factor of 2N f N e /N γ that arose because that front-end transformation related a diagram with an initial electron to one with an initial photon, and the 2N f N e /N γ reflected the different factors associated with averaging over initial flavors and helicities. In our case, the initial particle is always a gluon, so no such adjustment is necessary. Also, eqs. (H.13) and (H.14) of ref. [15] do not have the overall minus sign of our (2.7) above because they included a back-end transformation in addition to the front-end transformation. Note that those equations have also implemented x ↔ y in addition to the front-end transformation (2.7) above. The overall factor of 1 2 is included because of the loop symmetry factor associated with the (red) loop in theȳxȳx virtual diagram.
The only two virtual diagrams not covered by figures 9 and 10 are xyyx and xȳȳx. But these diagrams are related to each other by combined front-end and back-end transformations, as depicted in fig. 11. That means that transformations have given us a short-cut for determining all virtual diagrams except for one, which we take to be xyyx. Fortunately, the xyyx diagram has the same form as the QED diagram of fig. 12 previously computed in ref. [15], and the QED result can be easily adapted to QCD. One just needs to include QCD group factors associated with splitting vertices; use QCD instead of QED Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions; correctly account for identical-particle symmetry factors; and use QCD rather than QED results for the complex frequencies and normal modes associated with theq approximation to the propagation of the high-energy particles through the medium. Details of the conversion are given in appendix D 4.
We give more detail on implementing the above methods in appendix D, and final results for unrenormalized diagrams are given in appendix A [with σ ren =0 and σ bare =1 in section A 3].
C. UV divergences, renormalization, and running of α s The virtual diagrams of figs. 4 and 5 contain UV-divergent loops in the amplitude or conjugate amplitude. It may seem surprising that most of them can be related via figs. 9 and 10 to real double splitting (g → ggg) diagrams that involve only tree-level diagrams in the amplitude and conjugate amplitude. This is possible because we are working with time-ordered diagrams: individual time-ordered interferences of tree-level diagrams are UVdivergent even though the sum of all the different time-orderings is not. See section 4.1 of ref. [15] for more discussion of this point. In any case, the original calculations [21][22][23] of the g→ggg diagrams of figs. 1 and 2 discussed the UV divergence of each diagram and showed that they indeed canceled.
The corresponding divergences of the virtual diagrams, however, will not cancel. Indeed, they must conspire to produce the known renormalization of α s . Ref. [15] demonstrated how this worked out for large-N f QED, but the diagrammatics of renormalization of the QCD coupling is a little more complicated. We will also encounter a well-known annoyance of Light Cone Perturbation Theory (LCPT): individual diagrams will contain mixed UV-IR divergences that only cancel when the diagrams are summed together. 10

UV and IR regulators
We use dimensional regularization in 4−ǫ spacetime dimensions for UV divergences. However, we use the letter d to refer to the number of transverse spatial dimensions (2.9) For infrared divergences, we introduce a hard lower cut-off (p + ) min on light-cone momentum components p + . Hard momentum cut-offs complicate gauge invariance, but this is a fairly standard procedure in LCPT, since LCPT is formulated specifically in light-cone gauge A + =0. Note that p + is invariant under any residual gauge transformation that preserves light-cone gauge. It would of course be nicer to use a more generally gauge invariant choice of infrared regulator, but that would lead to more complicated calculations. 11 We will write our IR cut-off on longitudinal momenta p + as (p + ) min = P + δ (2.10) where P + is the longitudinal momentum of the initial particle in the double-splitting process and δ is an arbitrarily tiny positive number. 12 For consistency of IR regularization of the theory, this constraint must be applied to all particles in the process. For instance, in a g→ggg process where P + splits into daughters with longitudinal momenta xP + , yP + , and zP + , we require that the longitudinal momentum fractions x, y, and z all exceed δ.
(This automatically guarantees that internal particle lines in g→ggg diagrams also have p + > P + δ.) In a virtual correction to g → gg where P + splits into xP + and (1−x)P + , we must have x and 1−x greater than δ, but we must also impose that the momentum fractions of internal virtual lines are greater than δ as well. We'll see explicit examples below. With 10 For an example from calculations that are tangentially related to ours, see Beuf [26,27] and Hänninen, Lappi, and Paatelainen [28,29] on next-to-leading-order deep inelastic scattering (NLO DIS). For a description of the similarities and differences of our problem and theirs, see appendix B of ref. [15]. For a very early result on obtaining the correct renormalization of the QCD coupling with LCPT in the context of vacuum diagrams, see ref. [30]. 11 In particular, one might imagine using dimensional regularization for the infrared as well as the ultraviolet.
Unfortunately, the dimensionally-regulated expansions in ǫ that we currently have available [15,23] for the types of diagrams we need all made use of the fact that dimensional regularization was only needed for the ultraviolet. 12 A technicality concerning orders of limits: One should take the UV regulator ǫ → 0 before taking the IR regulator δ → 0. Taking δ → 0 first would be equivalent to using dimensional regularization for the IR as well as the UV, which is currently problematic for the reason given in footnote 11. this notation, the annoying mixed UV-IR divergences of LCPT are proportional to ǫ −1 ln δ, which is the product of a logarithmic UV divergence ǫ −1 and a logarithmic IR divergence ln δ.

Results for UV (including mixed UV-IR) divergences
We can read off the results for 1/ǫ divergences from the complete results given in appendix A. However, we will take the opportunity to be a little more concrete here in the main text by stepping through the calculation for one of the diagrams, but focusing on just the UVdivergent (1/ǫ) terms. Then we'll put the diagrams together to see the cancellation of mixed UV-IR divergences and the appearance of the QCD beta function coefficient β 0 .
Consider the first NLO g→gg diagram (yxxy) in fig. 9, which shows that diagram related by back-end transformation to the g→ggg diagram xyȳx. The 1/ǫ piece of the latter can be taken from ref. [23] and is [see appendix B of the current paper for more detail] where and (α, β, γ) are functions of x and y that represent various combinations of the helicitydependent DGLAP splitting functions associated with the vertices in the diagram. 13 In this section we use ≈ to indicate that we are only keeping 1/ǫ terms. Back-end transforming the above expression and swapping x↔y, as indicated in fig. 9, gives the corresponding result for the virtual diagram yxxy: where we have taken 2 Re(· · · ) to include the conjugate diagram as well. Doing similar calculations for the other crossed Class I diagrams (the top line of fig. 9), by using g→ggg results for xȳyx and xȳxy from ref. [23] and then transforming as in fig. 9, 13 Details of the definition of (α, β, γ) in terms of DGLAP splitting functions are given in sections 4.5 and 4.6 of ref. [21]. In order to make those definitions work with front-end transformations, one must additionally include absolute value signs as discussed after eq. (A23) of the current paper. (2.14) Since we are focused here just on the 1/ǫ pieces above, the integral may be done using the explicit d=2 expressions (A23) for (α, β, γ). But the combination (α + γ)z + (β + γ)(1−x)(1−y) appearing in (2.14) turns out to be dimension-independent in any case! (See appendix C.) Remember that for the crossed virtual diagrams, like all the Class I diagrams of fig. 4, taking x → 1−x generates other distinct diagrams that need to be included as well. So, do the y integral in (2.14), combine the result with x → 1−x [as in (1.5b) or (1.9)], and take the small-δ limit. This gives and P (x) is the DGLAP g→gg splitting function. A non-trivial feature of (2.15) is that the y integration in (2.14), combined with the addition of x → 1−x, gave a result proportional to the P (x) in (2.16). This is what will later make possible the absorption of 1/ǫ divergences by renormalizing the α s in the leading-order result. For the time being, however, note the unwanted mixed UV-IR divergence ǫ −1 ln δ in (2.15). Now turn to the sequential virtual diagrams. The sum 2 Re[xyxȳ + xxyȳ + xxȳy] of non-virtual sequential g→ggg diagrams shown in fig. 2 (together with their conjugates) represents the sum of all time orderings of a tree-level process and so does not give any net 1/ǫ divergence. 15 So there will also be no divergence in its back-end transformation, which fig. 9 shows is equivalent to the sum 2 Re[xyxy + xxyy + xxȳȳ] of three Class I sequential virtual diagrams. Nor will there be any divergence to its front-end transformation followed by the swap x ↔ y, corresponding by fig. 10 to the sum 2 Re[ȳxȳx +ȳȳxx + yyxx] of three Class II sequential diagrams. So none of these groups of diagrams generate a divergence.
What remains of figs. 4 and 5 is the Class I virtual diagram xyyx and the Class II virtual diagram xȳȳx, which are related to each other via fig. 11. As mentioned earlier, the result for 2 Re[xyyx] can be converted from the known result [15] for the similar QED diagram of fig. 12. The UV-divergent 1/ǫ piece of that QED result was 16 (2.17) The translation from a QED diagram to a QCD diagram is explained in our appendix D 4 and gives Our IR cut-off δ must now be included with the integration limits because, unlike QED, LPM splitting rates are (non-integrably) infrared divergent in QCD. The sgn M factors are included above because, even though M −1,x,1−x is positive for the xyyx diagram, this more general form is consistent with the front-end transformation we are about to perform. Since xyyx above is a Class I diagram, we need to also add in the other diagram that is generated by x → 1−x. Finally, the transformation of fig. 11 gives the remaining (Class II) diagram xȳȳx. 17 The sum of all three is Adding (2.15) and (2.19) gives the total UV divergence from virtual corrections to single splitting: y e ≡ y e /(1 − x e ). There was an overall sign error in eq. (F.42) of the original published version of ref. [15], which is treated correctly in the version above. 17 As discussed after eq. (A5), one must include an absolute value sign in the definition of P (x) in order to make it work with front-end transformations using our conventions.
The β 0 above is the same coefficient that appears in the one-loop beta function for α s = g 2 /4π: 22) where N f is the number of quark flavors. The N f term does not appear in (2.21) because we have not included quarks in our calculations, consistent with our choice to work in the large-N c limit (for N f fixed).
Note that the UV-IR mixed divergences have canceled between (2.15) and (2.19), as well as the ln x(1 − x) terms. These cancellations had to occur in order for the total divergence of the virtual diagrams to be absorbed by usual QCD renormalization, as we'll now see. 18

Renormalization
Following ref. [15], 19 we find it simplest to implement renormalization in this calculation by imagining that all diagrams have been calculated using the bare (unrenormalized) coupling and then rewriting (α s ) bare in terms of (α s ) ren . For the MS-renormalization scheme, that's α bare  18 There is something sloppy one might have tried in the preceding calculations that would have failed to produce the correct UV divergences, which we mention here as a caution to others because we unthinkingly tried it on our first attempt at this calculation. Suppose that we had set δ to zero in all the integration limits so that each IR-divergent integral we've done was divergent and ill-defined. Then suppose that in each integral we scaled the integration variable y so that each integral was now from 0 to 1, e.g.
1 0 dy f (1−x)y and similarly for x → 1−x. Now that the integration limits are the same, one could add together all the integrands for all the diagrams. The combined integral would be convergent but does not give the correct result (2.20). That's because one can get any incorrect answer by manipulating sums of ill-defined integrals. To properly regularize a theory, one must first independently define the cut-off on the theory (in this case the IR cutoff on longitudinal momenta) and only then add up all diagrams calculated with that cut-off. 19 Note that, because it is multiplied by 2/ǫ, we will need to use a d=2−ǫ formula for [dΓ/dx] LO in the last term above, as indicated by the subscript. We can now use (2.25) to regroup terms in (2.24) to write the LO+NLO g→gg rate in terms of MS renormalized quantities as One can see from (2.20) that the 1/ǫ poles indeed cancel in this renormalized [∆dΓ/dx] NLO . There are many equivalent ways to introduce the MS renormalization scale into the renormalization procedure outlined above. Following ref. [15], 21 we will introduce it by writing the dimensionful bare g 2 /4π in 4−ǫ spacetime dimensions as µ ǫ α bare s , where α bare s is the usual dimensionless coupling for 4 spacetime dimensions. As a result, every power of α s in our unrenormalized calculations comes with a power of µ ǫ which, if multiplied by a 1/ǫ UV divergence and expanded in ǫ, will generate the correct logarithms ln µ of the renormalization scale in our results, as we detail next.

Organization of Renormalized Results
Formulas for the NLO g→gg rate are given in appendix A 3. Because of the fact that multiple diagrams contribute to cancellation of 1/ǫ poles in ways that are not particularly simple diagram by diagram, we have organized our renormalized result for [dΓ/dx] NLO,ren g→gg slightly differently than the QED case of ref. [15], in a way that we will explain here.
Also, we would like to write renormalized formulas in appendix A 3 in a way that makes transparent the dependence on explicit renormalization scale logarithms ln µ. The running (2.22) of α s , plus the fact that the leading-order rate is proportional to α s , implies that the renormalized NLO rate must have explicit µ dependence in order to cancel the implicit µ dependence dα s /d(ln µ) = β 0 α 2 s of α s (µ) from the LO rate. In contrast, the NLO bare rate [∆ dΓ/dx] NLO,bare g→gg is proportional to (µ ǫ α s ) 2 , and so its divergence (2.20) generates (2.29) 21 See in particular the discussion of eq. (F.31) of ref. [15].
The difference between the ln µ terms of (2.28) and (2.29) is made up by the last term of the renormalization (2.27), as we'll now make explicit while also keeping track of all O(ǫ 0 ) pieces of the conversion.
To start, we need the d=2−ǫ dimensional result for the leading-order single splitting process, which appears in (2.27). We'll find it convenient to write this as 2π Here B(x, y) ≡ Γ(x) Γ(y)/Γ(x+y) is the Euler Beta function; we use the short-hand notations Ω 0 and M 0 for and the DGLAP g → gg splitting function P (x), given by (A5), is independent of dimension (see appendix C). Using (2.30), we rewrite the renormalized rate (2.27) as The first term [dΓ/dx] ren log of (2.33) contains the correct explicit ln µ dependence of (2.28 This is the meaning behind the notation we used back in (1.5). The notation is convenient because, for our final renormalized g→gg results listed in appendix A, the notation distinguishes the parts [∆ dΓ/dx] NLO of our results that are expressed in terms of y integrals, 23 like in (1.6), from the parts [dΓ/dx] LO above that are not.

III. IR DIVERGENCES IN ENERGY LOSS CALCULATIONS
We now discuss in detail how the IR behavior of various measures of the development of in-medium high-energy QCD parton showers depends only on the combination v(x, y) of virtual and real diagrams introduced in (1.11), for which power-law IR divergences cancel. In this section, ≈ indicates an equality that is valid at leading-log order.

A. General shower evolution
We start by looking generally at the evolution of the distribution of partons in such a shower. This will generalize, to NLO, similar methods that have been applied by Blaizot et al. at leading order [31,32]. 24 In what follows, let E 0 be the energy of the initial parton that starts the entire shower. We will let ζE 0 refer to the energy of some parton in the shower as the shower develops, and we will refer to the distribution of shower partons in ζ at time t as N(ζ, E 0 , t). Formally, the total number of partons remaining in the shower at time t is then 1 0 dζ N(ζ, E 0 , t), but this particular integral is IR divergent, not least because some fraction of the energy of the shower will have come to a stop in the medium (ζ=0) and thermalized by time t. However, one may also use N(ζ, E 0 , t) to calculate IR-safe characteristics of the shower, including N(ζ, E 0 , t) itself for fixed ζ > 0. 25

Basic Evolution Equation
The basic evolution equation to start with is (see appendix B for some more detail) 26 refers to the net rate (1.5) to produce one daughter of energy xE (plus any other daughters) via single splitting or overlapping double splitting from a parton of energy E. The total splitting rate Γ in the loss term is where the 1/2! and 1/3! are the final-state identical particle factors for g → gg and g → ggg.
The first and second terms in (3.2) are respectively loss and gain terms for N(ζ, E 0 , t).
The gain term corresponds to the rate for any higher-energy particle in the shower (energy ζE 0 /x) to split and produce a daughter whose energy is ζE 0 . To keep formulas simple here and throughout this discussion, we will not explicitly write the IR cut-off δ in integration limits.
By comparing (3.4) to (1.5), note that because of the different combinatoric factors involved in how [dΓ/dx] net is defined. This is related to the fact that (3.2) should not conserve the total number of partons: each g → gg should add a parton, and each g → ggg should add two partons. 27 The various pieces that go into the calculation of the right-hand side of the evolution equation (3.2) have various power-law IR divergences which cancel in the combination of all the terms. We now focus on identifying those divergences and showing how to reorganize (3.2) into an equivalent form where power-law IR divergences are eliminated from the integrals that must be done.

x → 0 or 1 divergences at leading order
To start, let's ignore NLO corrections for a moment and look at the leading-order version of (3.2): The leading-order rate [dΓ/dx] LO diverges as [see eq. (2.16) with ǫ=0]. Up to logarithmic factors, this divergence is the same for [dΓ/dx] LO (2.37) as well. This means that the integral (3.7) that gives the total rate Γ LO generates power-law IR divergences from both the x → 0 and x → 1 parts of the integration region. In contrast, the integral for the gain term in (3.6) runs from ζ>0 to 1 and so only generates a divergence from the x → 1 behavior. That means that we cannot get rid of the IR divergences simply by directly combining the integrands. However, if we first use the identical final- 27 One way to see this clearly is to over-simplify the problem by pretending that splitting rates did not depend on energy E, then integrate both sides of (3.2) over ζ, and rewrite then we can combine the loss and gain terms in (3.6) into Similar to (1.7), we have implemented the actual limits of integration here using step functions θ(· · · ) so that we may combine the integrands. Because of the θ functions, the integrand has no support for x → 0 and so no divergence associated with x → 0. Because we have combined the integrands, however, one can see that the integrand behaves like 1/(1−x) 1/2 instead of 1/(1−x) 3/2 (3.8) as x → 1 because of cancellation in that limit between the loss and gain contributions. So the form (3.10) has the advantage that the integral is completely convergent, and there are no IR divergences in this equation for any given ζ > 0.

y → 0 divergences at NLO
As discussed in section I B, g→ggg and NLO g→gg processes generate power-law IR divergences as the energy of the softest real or virtual gluon (whose longitudinal momentum fraction we often arrange to correspond to the letter y) goes to zero. We have already discussed how those power-law IR divergences cancel in the combination [∆ dΓ/dx] NLO net (1.7), which is the combination that appears in the NLO contribution to the gain term in the evolution equation (3.2). But the loss term involves a different combination Γ (3.4) of real and virtual diagrams, and so we must check that a similar cancellation occurs there.
Recalling that our NLO g → gg diagrams consist of our Class I diagrams ( fig. 4), their x → 1−x cousins, and our class II diagrams ( fig. 5), the NLO contribution to the total rate (3.4) is, in more detail, [Compare and contrast to (1.5b).] Fig. 13 shows the various integration regions corresponding to the different terms above and the limits of integration producing IR divergences (which is all of them).
We will now align the location of the IR divergences so that we can eventually combine the different integrals and eliminate power-law divergences. First, note by change x → 1−x of integration variables, the "(x → 1−x)" term in (3.11) gives the same result as the "virt I" term. Second, simultaneously use the x → 1−x and y → 1−y symmetries of Class II diagrams to divide the integration region of fig. 13c in half diagonally, giving For the NLO g→gg contributions, we now divide the integration region into (i) 0 < y < (1−x)/2 and (ii) (1−x)/2 < y < 1−x and change integration variables y → z = 1−x−y in the latter, similar to the manipulations used earlier to obtain (1.7). For the g→ggg contributions, note that permutation symmetry for the three final daughters (x, y, z) implies the integral over each of the six regions shown in fig. 14 is the same. We can therefore replace the integral over all six regions by three times the integral over the bottom two, depicted by the shaded region of fig. 15d. [We will see later the advantage of integrating over these two regions instead of reducing the integral to just one region.] Eq. (3.12) can then be written as with v and r defined as in (1.8). We will find it convenient to change integration variable x → 1−x in the first term and rewrite the equation as 14) The integration regions corresponding to the two terms in (3.14) are shown in fig. 15, where the only IR divergences correspond to y→0 or x→1.
The rationale for the last change was to convert x→0 divergences into x→1 divergences (the blue line in fig. 15), which we will later see then cancel similar x→1 divergences in the gain term of the evolution equation. For the moment, however, we focus only on the y→0 divergences of (3.14), depicted by the red lines in fig. 15. In the limit y → 0 (for fixed x), the integrand in (3.14) approaches where the first equality follows because g → ggg is symmetric under permutations of (x, y, z).
The right-hand side of (3.15) is the same combination as (1.11) but with x → 1−x. In fig.  8, we verified numerically that y −3/2 divergences (which generate power-law IR divergences when integrated) indeed cancel in this combination, leaving behind the double-log divergence 14: Equivalent integration regions for g → ggg corresponding to permutations of the daughters (x, y, z). The common vertex of these regions is at (x, y, z) = ( 1 3 , 1 3 , 1 3 ). shown in (1.11) [which happens to be symmetric under x → 1−x]. Interested readers can find non-numerical information on how the y −3/2 divergences cancel in appendix E.
One can now see why we did not replace the integral of r(x, y) over the two sub-regions shown in fig. 15 by, for example, twice the integral of just the left-hand sub-region (x < y < z). If we had done the latter, there would be no r term for x > 1/2 and so nothing would cancel the y −3/2 divergence of v(1−x, y) for x > 1/2. We had to be careful how we organized things to achieve our goal that the y integral in (3.14) not generate a power-law IR divergence for any value of x.
Next, we turn to our final goal for this section of showing that the integrals in the evolution equation for N(ζ, E 0 , t) can be arranged to directly avoid power-law IR divergences for the entire integration over both x and y.
4. x → 0 or 1 divergences at NLO By using (1.7), (3.10), and (3.14) in the shower evolution equation (3.2), we can now combine integrals to avoid all power-law divergences: (3.16c) We've previously seen that the LO piece S LO is free of divergences. And we've seen that the loss and gain terms of the NLO piece S NLO are each free of power-law divergences associated with y → 0 (with fixed x). Now consider divergences of S NLO associated with the behavior of x. The integrand in (3.16c) has no support as x → 0 (fixed y). And for x → 1 (fixed y), there is a cancellation between the loss and gain terms. So there is no divergence of S NLO associated with x → 0 or x → 1. 28 In summary, the only IR divergences coming from S NLO are the uncanceled double-log divergences associated with y → 0.

B. Absorbing double logs intoq and comparison with known results
Refs. [7][8][9] have previously performed leading-log calculations of overlap corrections and shown that the double-log IR divergences can be absorbed into the medium parameterq. We will now verify that the double-log piece of our results produces the same modification [36] ofq.

Double-log correction for [dΓ/dx] net
Let's start with the relatively simple situation of the [dΓ/dx] net introduced in section I B. From the discussion of (1.7) through (1.11), the double-log divergence of the NLO contribution to [dΓ/dx] net is given by 29 where we have re-introduced our sharp IR cut-off δ. (2.12) and (2.16)], the double-log correction above can be absorbed at this order by replacingq bŷ The corresponding leading-log modification ofq from earlier literature [7][8][9]36] is usually expressed in the final formq where L is the thickness of the medium and τ 0 is taken to be of order the mean free path for elastic scattering in the medium. In order to compare (3.19) and (3.20), we need to translate. First, for simplicity, we have been working in the infinite-medium approximation, which assumes that the size of the medium is large compared to all relevant formation lengths. Eq. (3.20) instead focuses on the phenomenologically often-relevant case where the width L of the medium is the formation time t form (x) associated with the harder splitting x. One may convert at leading-log level by considering the boundary case where Parametric substitutions like this inside the arguments of logarithms are adequate for a leading-log analysis. What remains is to translate between the use of two different types of cut-offs in (3.19) and (3.20): δ and τ 0 . To understand the effect of the cut-offs, it is useful to review where double logs come from in theq approximation, at first ignoring the cut-offs altogether. Parametrically, the IR double log arises from an integral of the form over the integration region shown in fig. 16a, given by 30 Using t form (y) ∼ yE/q for small y, these inequalities can be equivalently expressed as a range on y:q Now consider two different ways to evaluate the double logarithm (3.22). The first method is to add a lower cut-off τ 0 on ∆t, as in fig. 16b. Using (3.23b), that's Alternatively, adding a lower cut-off δ on y as in fig. 16c, using (3.21), and assuming When we extract just the double log dependence ln 2 δ on the parameter δ, there is no difference (for fixed x) at leading-log order between ln 2 (δ/x) and ln 2 δ. At that level, comparison of (3.24) and (3.25) gives the leading-log translation 30 Using (3.21) and t form (ξ) ∼ ξE/q for small ξ, (3.23a) can be put in the form y E/xq ≪ ∆t ≪ yE/q presented in eq. (9.3) of ref. [21] for y ≪ x ≤ z. The equivalence, in turn, with notation used in some of the original work on double logs in the NLO LPM effect is discussed in appendix F.1 of ref. [21].
between IR-regularization with τ 0 and δ. Applied to the standard double log result (3.20), this translation exactly reproduces the double log behavior (3.19) of our own results. We will return to the x dependence of (3.25) when we later examine sub-leading single-log corrections in section IV.
Our δ is simply a formal IR regulator. In contrast, there is a plausible physical reason for using the elastic mean free path τ 0 as an IR regulator at the double log level: Theq approximation used throughout our discussion and earlier literature is a multiple-scattering approximation that requires long time periods compared to the mean free time between collisions. However, beyond leading-log order, the use of a τ 0 cut-off would be problematic for full NLO calculations. In our calculations, a τ 0 cut-off would interfere with the correct UV-renormalization of α s , which comes from ∆t → 0 (and small enough time scales that evenq-approximation propagators faithfully reproduce vacuum propagators). So in this paper we have just chosen the formal IR regulator, δ, that seemed most convenient for our calculations.
In order to use IR-regulated results for NLO splitting rates, one must either compute quantities that are IR-safe in theq approximation or else make an appropriate matching calculation for soft emission that takes into account how the QCD LPM effect turns off for formation lengths τ 0 .

Physics scales: What if you wanted to take δ more seriously?
Though we are simply taking δ as a formal IR cut-off for calculations involving theq approximation, we should mention what the physics scales are where ourq-based analysis would break down if one used our results for calculations that were sensitive to the value of δ. The situation is complicated because there are potentially two scales to consider, indicated in fig. 17. We have given parametric formulas for those scales for the case of a weakly-coupled quark-gluon plasma. One may translate to a strongly-coupled quark-gluon plasma, in both the figure and the discussion below, simply by erasing the factors of g.
Parametrically, the mean free time between (small-angle) elastic collisions with the medium is τ 0 ∼ 1/g 2 T , andq is ∼ g 4 T 3 . Using the limits (3.23b) on y, as well as (3.21) and t form(x) ∼ xE/q, one then finds for ∆t ∼ τ 0 the corresponding soft gluon energies yE indicated in the figure.
Our formalism breaks down for yE smaller than the lower limit yE ∼ T because gluons of energy T cannot be treated as high-energy compared to the plasma. Note that if one correspondingly chose δ ∼ T /E without also constraining ∆t, then the resulting double log region would be larger than has been conventionally assumed in the literature. In contrast, if one chose δ ∼ xT /E, corresponding to the other red line in fig. 17, then one would guarantee that ∆t τ 0 but the resulting double log region would be smaller than the one used in the literature. There is no choice of δ alone that corresponds to the traditional shaded region of fig. 17.

Double-log correction for shower evolution equation
The gain term of the shower evolution equation (3.2) depends only on the combination [dΓ/dx] net of rates, and so the same redefinition (3.19) will absorb the double logarithmic divergence. One expects that this must also work for the loss term in (3.2), which depends on the combination Γ, but we should make sure. Since we found that only y → 0 ultimately contributes to the double logarithm in our later version (3.16) of the evolution equation, we can focus on the y→0 behavior of the NLO loss term for fixed x, which corresponds to the y→0 behavior of the integrand of (3.14) for ∆Γ NLO . Using (3.15) and (1.11), the double log generated by the y integration in (3.14) is When combined with the leading-order rate Γ LO given by (3.7), we have In the literature, it is common to discuss energy loss per unit length (dE/dL) for a highenergy particle. This makes sense only if one can unambiguously identify the original particle after a process that has degraded its energy. For many applications of the LPM effect, the energy loss occurs by radiation that is soft compared to the initial particle energy E, and so one can identify the particle afterwards as the only one that still has very high energy. In this paper, however, we have been focused on the case of a very thick medium (thick compared to formation lengths). In that case, hard bremsstrahlung is an important aspect of energy loss. If the two daughters of a splitting have comparable energies, it becomes more difficult to say which is the successor of the original. For a double-splitting process beginning with a quark, one can unambiguously (for large N c ) choose to follow the original quark. But, for processes that begin with g→gg, the distinction is less clear.
One possibility might be to formally define dE/dL for g→ gg processes by always following after each splitting the daughter gluon that has the highest energy of the two daughters. Unfortunately, this procedure is ill-defined when analyzing the effect of overlapping formation times on successive splittings. Consider the interference shown in fig. 18 of two different amplitudes for double splitting g → gg → ggg. For each amplitude, the red gluon line shows which gluon we would follow by choosing the highest-energy daughter of each individual g→gg splitting. The two amplitudes do not agree on which of the final three gluons is the successor of the original gluon. That's not a problem if the individual splittings are well enough separated that the interference can be ignored, i.e. if formation lengths for the individual splittings do not overlap. But since we are interested specifically in calculating such interference, we have no natural way of defining which gluon to follow. This is why we have avoided dE/dL and focused on more general measures of shower evolution.
The above argument generalizes to g → ggg points made in ref. [10] about e → γe →ēee, q → gq →qqq and q → gq → ggq. However, in those cases, ref. [10] noted that dE/dL was nonetheless well-defined in the large N f or N c limits respectively. In contrast, the g→ggg interference shown in fig. 18 is unsuppressed in the large-N c limit.

D. Similar power-law IR cancellations
LPM splitting rates and overlap corrections scale with energy like q/E, up to logarithms. For situations where rates are proportional to a power E −ν of energy, ref. [10] discusses how to derive relatively simple formulas for the stopping distance of a shower, and more generally formulas for various moments of the distribution of where the energy of the shower is deposited. Those formulas can also be adapted to the case where the rates also have single-logarithmic dependence E −ν ln E. This is adequate for analyzing stopping distances for QED showers [10], but the application to QCD, which has double logs, is unclear. But even for QCD, one can use those stopping length formulas as yet another context in which to explore the cancellation of power-law IR divergences. See appendix F for that analysis.

A. Numerics
In (1.11) and section III B 1, we extracted the known IR double logarithm from the slope of a straight-line fit to the small-y behavior of our full numerical results when plotted as v(x, y) + 1 2 r(x, y) vs. ln y, as in fig. 8. The sub-leading single-log behavior can be similarly found, for each value of x, from the intercept of that straight-line fit. Specifically, refine (1.11) to include single-log effects by writing v(x, y) Here, the y −1 ln y term generates the known double-log behavior ∝ ln 2 δ after integration over y, and the new s(x) y −1 term allows for additional single-log behavior ∝ ln δ. Then the combination (4.1) behaves at small y like v(x, y) + 1 2 r(x, y) The right-hand side represents the straight line fit of fig. 8, and the intercept of that fit at ln y = 0 gives −s(x). Our numerical results for s(x) are shown by circles in fig. 19. Note that s(x) is not symmetric under x → 1−x. That's because we defined v(x, y) in (1.8a) to contain Class I virtual diagrams but not their x → 1−x cousins. We do not have anything interesting to say about the precise shape of s(x) itself. But we can get to something interesting if we note that our original discussion (1.11) of the smally behavior of v(x, y) + 1 2 r(x, y) was in the context of [dΓ/dx] net , where v(x, y) of (1.10). For this combination, the single log piece corresponds to twice the averagē Thiss(x) is depicted by the diamonds in fig. 19. And even though we currently have only numerical results fors(x), we will be able to make some interesting observations about its form by comparing our numerics to an educated guess that we will discuss in a moment.
[dΓ/dx] net , and thuss(x), also appears in our other discussions of IR behavior, such as the gain term in the evolution equation (3.2) for the gluon distribution N(ζ, E 0 , t). The loss term of that equation depends on the total rate Γ, which treats the two identical daughters of g → gg processes x and 1−x on an equal footing. 31 Sos(x) is the relevant function for single log divergences, regardless of the fact that we found it convenient to rewrite Γ in (3.16) in a way that obscured the x ↔ 1−x symmetry of g→gg so that we could make more explicit the cancellation of power-law IR divergences. 32

B. Educated guess for form ofs(x)
Let's now return to the issue of x dependence in the translation of the standard double log result ln 2 (L/τ 0 ) in (3.24) to the ln 2 δ of our calculations in (3.25). Previously, when we compared the two, we ignored the x dependence of the ln 2 (δ/x) in (3.25). Now keeping 31 As was true for [dΓ/dx dy] net , the r(x, y) contribution representing g → ggg is symmetric in x ↔ z ≡ 1−x−y rather than x ↔ 1−x, but the difference is unimportant in the y→0 limit we are using to extract IR divergences. More specifically, the difference between r(x, y) = r(1−x−y, y) and r(1−x, y) is parametrically smaller as y→0 than the 1/y terms responsible for the single-log IR divergence under discussion. 32 If desired, one could achieve both goals by replacing the integrand in (3.16) by its average over x ↔ 1−x.
track of that x dependence, the translation (3.26) becomes Here we assume x < 1−x, and the arguments of the double logarithms are only parametric estimates. Rewrite the right-hand side of (4.6) as ln 2 ∆ with ∆ ∼ δ/x. For x ≪ 1, this parametric relation suggests that ∆ ≃ #δ/x for some proportionality constant #. So (4.6) suggests that a more precise substitution for x ≪ 1 would be Eq. (4.7) contains information about the small-x dependence of the coefficient of the subleading, single IR-logarithm ln δ.
In a moment, we will attempt to generalize to a guess of the behavior for all values of x, but first let's see how (4.6) compares to our numerics. Consider the logarithms arising from a symmetrizeds version of (4.2), whose integral over y would be proportional to − δ dy ln y +s(x) y = 1 2 ln 2 δ +s(x) ln δ + (IR convergent). (4.8) Comparison of (4.7) with (4.8) suggests that where c = ln # is a constant that is not determined by this argument and must be fit to our numerics. The dashed blue curve in fig. 19 shows (4.9) with * c = 9.0 (4.10) on the graph of our full numerical results. The form (4.9) works well for small x.
To make an educated guess for the full x dependence ofs(x), we need to replace (4.9) by something symmetric in x ↔ 1−x. The formation time t form (x), related to the harmonic oscillator frequency Ω 0 of (2.31) by is symmetric in x ↔ 1−x and plays a major role in the LPM effect. So, even though our arguments about double logs have only been parametric, let us see what happens if we guess that the 1/x in (4.9) is arising from the small x behavior of (4.11), and so we replace (4.9) bys (x) = ln −1 + 1 x C. How well does the educated guess work?
As the figure shows, (4.12) captures the x dependence of the single log coefficients(x) very well. However, it is not quite perfect. To see the discrepancies, one may use (4.2) together with (4.12) to extract from our numerical results for v(x, y) + 1 2 r(x, y) the best choice c(x) of c for each individual value of x: (4.13) If the guess (4.12) for the form ofs(x) were exactly right, then c(x) would be an xindependent constant. But fig. 20 shows a small variation of our c(x) with x. Our educated guess is a good approximation but appears not to be the entire story for understanding IR single logs. The variation of c(x) in fig. 20 is the reason that we have not bothered to determine the small-x value of c in (4.9) to better precision than (4.10).
We should note that the value of c will be IR-regularization scheme dependent. If we had regulated the IR with a smooth cut-off at p + ∼ P + δ instead of a hard cut off, a different value of c would be needed to keep the physics the same on the right-hand side of (4.8) with the different meaning of δ.

V. THEORIST ERROR
The results presented in Appendix A for overlap effects on double splitting calculations represent the culmination of a very long series of calculations [15,[21][22][23] that required addressing many subtle technical issues as well as many involved arguments computing expansions in ǫ for novel dimensionally-regulated quantities. In the absence of calculations by an independent group using independent methods, a natural worry must be whether somewhere our group might have made a mistake that would noticeably affect our final results. We refer to this possibility as "theorist error," in contrast to "theoretical error" estimates of uncertainty arising from the approximations used.
Though we cannot absolutely guarantee the absence of theorist error, we think it useful to list a number of cross-checks and features of our calculations. Some of these check our treatment of technical subtleties of the calculation.
1. The power-law IR divergences computed for real and virtual diagrams in theq approximation cancel each other, as discussed in this paper. Sub-leading IR divergences, which do not cancel, correctly reproduce the IR double log [36] known from previous, independent calculations [7][8][9] that analyzed overlap effects in leading-log approximation.
2. Our calculation generates the correct 1/ǫ UV divergences for the known renormalization of α s . This includes the cancellation of mixed UV-IR divergences, which is one of the subtleties of Light-Cone Perturbation Theory.
3. In the soft limit y ≪ x ≪ 1 of g → ggg, crossed [21] and sequential [22] diagrams give contributions to ∆Γ/dx dy that behave like ln(x/y)/xy 3/2 . But the logarithmic enhancement of these 1/xy 3/2 contributions cancels when all g→ggg processes are added together, reassuringly consistent with the Gunion-Bertsch picture presented in appendix B of ref. [22]. When our formalism is applied instead to large-N f QED [15], the analogous logarithm does not cancel. In that case, its coefficient reassuringly matches what one would expect from DGLAP-like arguments, as explained in section 2.2.3 of ref. [15]. 4. One of the technical subtleties of our methods has to do with identifying the correct branch to take for logarithms ln C of complex or negative numbers, which may arise in dimensional regularization, for example, from the expansion of a C ǫ . See section 4.6 and appendix H of ref. [23], as well as appendix H.1 of ref. [15], for examples where the determination of the appropriate branch requires care. Making a mistake of ±2πi in the evaluation of a logarithm would generally have a significant effect on our results. But we do have some consistency checks on such "π terms" that result from the logarithm of the phases of complex numbers in our calculation. One check is illustrated by appendix E, where π terms associated with individual diagrams must all cancel as one part of the cancellation of IR power-law divergences. A different, somewhat indirect cancellation test of π terms generated by dimensional regularization is given in appendix D of ref. [23].

5.
Here is another test of an O(ǫ 0 ) term in the expansion of dimensional regularization of a UV-divergent diagram. Recall that both g→ggg and NLO g→gg processes have power-law IR divergences of the form δ dy/y 3/2 ∼ δ −1/2 , where the power law y −3/2 matches a physical argument given in section I.D of ref. [22]. In the calculation of divergent diagrams, the UV-sensitive piece of the calculation is isolated into what are called "pole" pieces in refs. [15,[21][22][23] and in appendix A. These pole pieces are evaluated analytically with dimensional regularization and yield 1/ǫ divergences plus finite O(ǫ 0 ) contributions. The remaining UVinsensitive contributions to the diagrams are evaluated with numerical integration. For some of the crossed virtual diagrams (top line of fig. 4), both the O(ǫ 0 ) pole piece and the UV-insensitive numerical integral 33 turn out to have spurious IR divergences that are more IR divergent than the power-law divergences we have discussed. However, they also turn out to exactly cancel each other. For example, in appendix E 4, we show how the integral associated with 2 Re(xyȳx) has an unwanted dy/y 2 ∼ δ −1 divergence from y→0 that is canceled by the O(ǫ 0 ) piece of the UV-divergent pole term. 34

VI. CONCLUSION
The results of this paper (combined with those of earlier papers) are the complete formulas in appendix A for the effects of overlapping formation times associated with the various g→ggg and g→gg processes of figs. 1-5. But there are still missing pieces we need before we can answer the qualitative question which motivates this work: Are overlap effects small enough that an in-medium shower can be treated as a collection of individual high-energy partons, assuming one first absorbs potentially large double logarithms into the effective value ofq?
First, for a complete calculation, we will also need processes involving longitudinal gluon exchange and direct 4-gluon vertices, such as in fig. 7. The methods for computing those diagrams are known, and so it should only take an investment of care and time to include them.
More importantly, our results as given are double-log IR divergent. The known double-log IR divergence can easily be subtracted away from our results and absorbed into the effective value ofq reviewed in section III B 1. However, this potentially leaves behind a sub-leading single-log IR divergence. We've seen from numerics that much of those single-log divergences can also be absorbed intoq eff by accounting for the x dependence of the natural choice of scale for the double-log contribution toq eff , but there remains a smaller part of the single-log IR divergences that is not yet understood. In order to make progress and understand the structure of the single logarithms, we hope in the future to extract analytic (as opposed to numerical) results for them from our full diagrammatic results. We have also not yet determined whether diagrams involving longitudinal gluon exchange, which have so far been left out, contribute to IR single logarithms.
It would be extremely helpful, both conceptually and as a check of our own work, if someone can figure out a way to directly and independently compute the sub-leading singlelog IR divergences without going through the entire complicated and drawn-out process that we have used to compute our full results. 33 In formulas, the pole piece of the crossed virtual diagrams corresponds to eq. (A58) for A pole virt Ic . whereas the UV-insensitive piece is the integral shown in (A55). For more details on exactly how the pole piece is defined, see appendix D. 34 This is unrelated (as far as we know) to a different class of cases, where individual diagrams have unwanted IR divergences that are only canceled by similar divergences of another diagram. See the two pairs of dz/z 5/2 divergences in Table I in appendix E.

Note added
After this paper was published, we found an error in eq. (A37) for A pole seq (x, y), which gave incorrect iπ terms when the arguments x or y are negative and so generated an incorrect result when front-end transformed for evaluation of type II virtual sequential diagrams. The correct version is derived in appendix A of ref. [39]. Here, we have chosen to leave (A37) and figs. 19 and 20 as in the original publication but will describe the changes. The correct formula for A pole seq (x, y) that works with front-end transformations is quoted in a footnote below (A37). In this appendix, we collect final results for the elements contributing to the leading-order g→gg rate, its NLO corrections, and the g→ggg rate: Throughout this appendix, we define as in the main text. We remind readers that in this paper we have not included diagrams involving 4-gluon vertices or instantaneous interactions via longitudinal gauge boson exchange, such as the examples of fig. 7. 1. Leading-order splitting rate a. d=2 transverse spatial dimensions In our notation, the leading-order g→gg rate is with and the g→g DGLAP splitting function Here and throughout this paper, our P (x) is just the function above and does not include the pieces of the usual DGLAP splitting function used to include the effect of virtual diagrams. In particular, the 1/(1−x) above is just the ordinary function 1/(1−x) and not the distribution 1/(1−x) + , and our P (x) above does not contain a δ-function term δ(1−x). When we need to deal with virtual diagrams in this paper, we will do so explicitly. The absolute value signs in (A5) may seem redundant since the absolute value is taken of a quantity that is manifestly positive for 0 < x < 1. They are included so that our definition of P (x) works with front-end transformations, for the same reasons described after (A23) below.
For the sake of later formulas for virtual corrections, it will be helpful to also express the above result in terms of the xx diagram of fig. 3 as where [dΓ/dx] ren log is given by (A50).

g → ggg rate
For the diagrams considered in this paper, we have where the first term represents the crossed diagrams of fig. 1 and the second term the sequential diagrams of fig. 2. A summary of the formulas for these rates appears in appendix A of ref. [20]. We will also present them here (i) for convenient reference in this paper, especially since many of the new formulas we need are related, (ii) because some minor modifications are needed [to (A23) and (A46) below] to make the formulas work in a simple way with front-end transformations, (iii) because we've rewritten some old formulas [such as (A15)] in a way that makes clearer their relation to some new formulas [such as (A58)], and (iv) to include some notational definitions [such as (A20) and (A44)] that were omitted from the summary in ref. [20].

a. Crossed Diagrams
Here we collect the result for the crossed diagrams [21] as corrected by ref. [23]. A brief summary of the interpretation of each piece below can be found in section VIII of ref. [21].
A(x, y) = A pole (x, y) + Note that the (α, β, γ) used in the definition (A12) of B are implicitly functions (A46) of the arguments x and y of B(x, y, ∆t) [with z ≡ 1−x−y]. This is important in formulas such as (A11), where in some terms those local arguments are replaced by other variables. Eq. (A23) gives (α, β, γ) for d=2. However, as explained in appendix C (which gives the more general formulas for d=2−ǫ), the d=2 formulas for (α, β, γ) are all that is needed here in appendix A.
The absolute value signs in (A23) may seem unnecessary since g → ggg processes have parton longitudinal momentum fractions x, y, z, 1−x, 1−y all positive. The advantage of including an absolute value sign around every such parton momentum fraction [which is equivalent to the use of absolute value signs in (A23)] is that they make front-end transformations like (2.7) work in a simple way, despite the fact that the front-end transformation replaces x by a negative number. 35 Theq → 0 limit for the vacuum piece in (A13) corresponds to taking all Ω's to zero and so making the replacements where 1 is the identity matrix. For numerical evaluation, one must take care that the above takes X yȳ → 0 and so

4-particle frequencies and normal modes
Here we collect formulas for the large-N c frequencies and normal modes associated with 4-particle propagation (section V.B of ref. [21]). 35 They are the QCD version of the absolute value signs used in eq. (A22) of ref. [15], which are discussed in footnote 38 of ref. [15]. One could alternatively dispense with the absolute value signs in the QCD case (A23) above by noting that negating x in that formula would, without absolute value signs, simply introduces a common overall minus sign in the values of (α,β,γ), which could be accounted for by modifying the sign of the front-end transformation formula (2.7). We've chosen to introduce the absolute value signs, however, so that our overall sign convention for front-end transformations will be the same as it was in the QED case of ref. [15].

c. Sequential Diagrams
Here we collect the result for the sequential diagrams [22]. A brief summary of the interpretation of each piece below can be found in section III of ref. [22]. Symbols such as Ω ± or a y , which are written in the exact same notation as symbols defined above, are given by their definitions above.

NLO g → gg rate
With regard to renormalization, we are going to make our summary formulas in this subsection do double duty by introducing a variable σ ren ≡ 1, for renormalized results; 0, for unrenormalized results, (A47) and its complement σ bare ≡ 1 − σ ren .
In the renormalized case (σ ren =1, σ bare =0), the α s in the leading-order splitting rate (A3) is MS-bar renormalized α s with renormalization scale µ, and we have chosen to group all of the µ-dependence of the NLO diagrams into the term shown explicitly in (A50) below, as discussed in section II C 4. In the unrenormalized case (σ ren =0, σ bare =1), the α s in the leading-order splitting rate (A3) is instead the bare α s , and we show the 1/ǫ and ln µ dependence of the NLO diagrams individually for each diagram.
Along the lines discussed in section II C 4, we write In what follows, we will further subdivide Class I diagrams into what we call (Class Ic) crossed virtual diagrams, given by the first row of fig. 9 plus conjugates; (Class Is) back-end sequential virtual diagrams, given by the remaining three diagrams of fig. 9 plus conjugates; and 2 Re(xyyx), given by the last diagram of fig. 4 Above, C and (α, β, γ) are the same as (A13) and (A23) for g→ggg crossed diagrams.
Note: The shorthand notation Ω 0 (A4) used above is the same as the Ω −1,x,1−x (A20) also appearing above, but we have used the latter to make explicit the similar structures of the three terms in (A58).

b. Sequential Virtual Diagrams
Above, A seq is the same as (A32) for sequential g→ggg diagrams. See appendix D 3 for alternative ways to write (A60) and for comments concerning the physical meaning of the ∆t integration variable of (A32) in the context of (A60).
are the same as the I seq n of (A39) except that the (X, Y, Z) seq there are replaced by See appendix D 5 for alternative ways to write (A69) and for comments concerning the physical meaning of the ∆t integration variable.
where ζ ′ E 0 > ζE 0 is the energy of a particle in the shower that decays into a daughter carrying fraction x of the parent's energy. The δ function requires that the daughter's energy xζ ′ E 0 match the energy ζE 0 we are looking for, and all possibilities for ζ ′ and x are integrated over. Using the δ function to do the ζ ′ integral gives the gain term in (3.2).

Eq. (E18):
The desired integral is convergent, but it will be useful to integrate the two terms separately. We must introduce a regulator to split up the integration because the integral of each term by itself is divergent. So consider the more general convergent integral and follow logic similar to dimensional regularization. By scaling arguments, similar to dimensional regularization, the integral of any power must be zero. For example, and then differentiating this result with respect to ǫ gives Writing ln(aτ ) = ln a + ln τ , (B4) and (B5) then give and so the first term in (B3) integrates to zero with this regularization. We are left with Consider ǫ > 1 (for which this integral is convergent), and then later analytically continue to ǫ = 0. We can rewrite (B7) as From eq. (3.527.1) of Gradshteyn and Ryzhik [37], 36 where ζ is the Riemann zeta function. Eq. (B8) then gives the desired result for our original integral by taking the limit ǫ=0. To calm any doubts about this derivation, one may simply check the answer numerically.
Eqs. (A23) and (A46) present d=2 results from refs. [21] and [22] for (α, β, γ) and (ᾱ,β,γ), which are various combinations of helicity-dependent DGLAP splitting functions that arise in calculations of g → ggg diagrams. However, in this paper, we use these same quantities in the calculation of virtual diagrams for g → gg, which are UV-divergent. So one might expect that when an α or β or γ is multiplied by a divergent 1/ǫ, then we need to know the O(ǫ) corrections to (α, β, γ) in order to calculate the finite pieces of our g→gg virtual diagrams, similar to what happens for QED in ref. [15]. 37 In this appendix, we present d=2−ǫ results for (α, β, γ) and (ᾱ,β,γ). However, we will see that, in the final results of appendix A, (α, β, γ) and (ᾱ,β,γ) only appear in combinations where the O(ǫ) pieces cancel, and so the original d=2 results are all that are actually needed there.
The first important fact is that the helicity-averaged g→gg DGLAP splitting function P (x) given in (A5) does not depend at all on dimension and so has no O(ǫ) correction. [See, for example, eq. (17) of ref. [38], which one may verify independently.] This lack of dependence on dimension is special to helicity-averaged g→gg splitting. Splittings involving quarks do depend on dimension, but we do not consider those in the large-N c limit of gluon-initiated showers considered in this paper.
For the particular combinations of helicity-dependent splitting functions that we need, we found it easiest to do the calculation from scratch. The helicity basis is unwieldy in general dimensions since there are no longer simply two helicities ±, and we find it simpler to do the calculation in a basis of linear polarizations. Other than that, we will follow the same notation and normalization conventions and derivations that were used for the d=2 case in sections 4.5 and 4.6 of ref. [21]. (See also appendix C of ref. [21].) Following ref. [21], we write our splitting vertex matrix elements in the form where and the p's represent transverse momentum, with p i = p j + p k . The T i→jk factor above implicitly depends on the polarization, longitudinal momentum fractions, and color states of the parent i and daughters j, k. The color factor is the T color i→jk above, which is −if abc for g → gg.
One may then extract the splitting functions P i→jk from the corresponding matrix elements in the nearly-collinear limit, and it's easiest to do this by temporarily choosing the axes so that the parent has transverse momentum zero: p i = 0 above. Then define q ≡ p j = −p k . One can calculate that the matrix element M for the three-gluon interaction is given by Here, capital roman letters I, J, K run over 1, 2, · · · , d=2−ǫ and index a basis for the (linear) transverse polarization states of the particles i, j, k. ξ j = x j /x i and ξ k = x k /x i are the longitudinal momentum fractions of the two daughters relative to their immediate parent in g → gg. The x's are longitudinal momentum fractions of the various particles in this one particular g→gg splitting relative to the original particle that initiated the entire double splitting process. In the nearly-collinear limit relevant to high-energy bremsstrahlung, the energies of the particles in this g→gg splitting are then where E is the energy of the original particle that initiated the double-splitting process. We bring this up in order to match conventions with the analysis in ref. [21]. That analysis used non-relativistic normalization of states, and so the desired matrix element is related to the more conventional M rel above by With our temporary convention that p i = 0, we have P jk = (x k +x j )q = x i q. Then comparison of (C1) with (C4) gives the components of P i→jk to be In this appendix, we will assume that all the (x i , x j , x k ) are positive and will not bother with the absolute value signs that were included in (A23) to be consistent with front-end transformations.
We can now use (C5) in the definition of the combinations (α, β, γ) in eqs. (4.37-38) of ref. [21], which is α(x, y) δn n δm m + β(x, y) δnmδ nm + γ(x, y) δn m δ nm [Here, we've indexed the possible linear polarization states using the letter I, whereas in ref. [21] the helicity basis was used, indicated by the letter h there.] Plugging (C5) into the right-hand side of (C6) and doing all the sums over polarization indices for d transverse dimensions, we can then extract from (C6) the results for (α, β, γ). For d = 2, the results are given in (A23) here and were originally presented in ref. [21]. For general d, we find The (d−2)/d terms above cancel in the combination which is the only combination that appears multiplying a UV-divergent 1/ǫ in our results summarized in appendix A [see (A58)]. For that reason, there is no problem with just using the d=2 values (A23) in appendix A. A similar procedure determines (ᾱ,β,γ), which are defined by eqs. (E.2,E.3) of ref. [22] asᾱ (x, y) δn n δm m +β(x, y) δnmδ nm +γ(x, y) δn m δ nm One can check that these results satisfy the QCD version of the identity of eq. (F32) of ref.
[15]: 38ᾱ remembering that for the case of g → gg, the polarization-averaged splitting functions P (d) (x) do not in fact depend on dimension d.
In our summary of results in appendix A, (ᾱ,β,γ) either appear in formulas where there are no UV-divergent 1/ǫ factors, or else only appear implicitly in d-independent combinations like P (x) P (· · · ) in (A66) and (A70). So the general-d formulas (C10) are not necessary for our results.

Appendix D: Details on transforming previous work to NLO g→gg diagrams
In this appendix, we give more detail about computing NLO g→gg diagrams. Since many of those diagrams are transformations of g→ggg diagrams, we start with the latter.

Prelude: g→ggg Crossed Diagrams
Though previous work [21,23] has calculated g→ggg processes with dimensional regularization, those calculations were complete only for sums of crossed diagrams for which UV divergences 1/ǫ canceled (as they must for tree-level processes). The transformations to virtual crossed diagrams in fig. 9 do not involve such UV-canceling collections of g→ggg diagrams, and so we now need complete results for individual g→ggg crossed diagrams. Consistently combining calculations of UV divergences with finite numerical integrals requires going slightly beyond what was done in ref. [23], and here we will organize the calculation using the methods developed in ref. [15].
In our calculations, UV divergences arise as ∆t→0 divergences of single integrals ∞ 0 d(∆t) F (∆t) of some function F (∆t). The full integrals are complicated enough that we do not know how to do them analytically. As explained in section 4.3.2 of ref. [15], our method for isolating the UV divergences and combining them with numerical integration is to rewrite where F d (∆t) is the integrand in dimensional regularization for d=2−ǫ transverse spatial dimensions. Above, D 2 (∆t) is any convenient function that • matches the divergence of F 2 (∆t) as ∆t → 0; • falls off fast enough as ∆t → ∞ so that ∞ a d(∆t) D 2 (∆t) will converge for non-zero a; • is simple enough that ∞ a d(∆t) D 2 (∆t) can be performed analytically. The last integral in (D1) is convergent and can be performed numerically. The first term can be found analytically by simplifying the otherwise complicated integrand F d (∆t) by expanding it in small ∆t. The scare quotes around the limit "a→0" in (D1) mean that a→0 should be taken after the ǫ→0 limit. The exact choice of D 2 does not matter: the total (D1) will be the same.
Let's focus on the xyȳx diagram in fig. 1. The d=2 integrand for 2 Re(xyȳx), corresponding to F 2 (∆t) above, can be taken from ref. [21] and corresponds to a piece of our eqs. (A11-A12): F 2 (∆t) = 2 Re C(−1, y, z, x, α, β, γ, ∆t), with C given by (A13-A14). The small ∆t behavior of this result is given by eq. (5.46) of ref. [21] as where Following similar choices made in ref. [15], we could now take D 2 (∆t) to be, for example, the right-hand side of (D3) with the replacements 39 which has the same small-∆t behavior but falls off faster as ∆t → ∞. However, for the presentation in this paper, it will be less cumbersome to just wait until we have assembled all the other elements for the sum of crossed virtual diagrams and then choose a single overall D 2 appropriate to that sum. The information about 2 Re(xyȳx) we will keep track of for now is (i) the d=2 formula (D2) for its integrand and (ii) the first integral in (D1), which integrates over small times. The latter, dimensionally-regulated integral is given by ref. [ (up to terms that vanish as ǫ→0), and As discussed in section 6 of ref. [21], the other crossed g→ggg diagrams can be obtained by various substitutions: where the changes to (x 1 ,x 2 ,x 3 ,x 4 ) are also applied to our formulas (A19) defining (Ω, M) i and (Ω, M) f .

Crossed Virtual Diagrams
We now obtain results for the crossed virtual diagrams from the preceding expressions by using the transformations of fig. 9.

a. (∆t < a)[1] terms
Let's first focus on the "(∆t < a) [1]" terms, which trace back to (D7), using (D10) when relevant. We find 2 Re dΓ dx dy ∆t<a [1] virt Ic where the contributions from individual diagrams (in some cases complex conjugated) are The 1/ǫ pieces of these formulas are the divergences (2.13) presented in the main text. The subscript "virt Ic" in (D11) stands for "virtual crossed diagrams" (which are all a type of Class I diagram), as in (A53). The yxyx is a little different from the other diagrams above because it is the only one that involves a front-end transformation. Fig. 9 shows that 2 Re(yxyx) is given by a front-end transformation of 2 Re(xȳxy) followed by x ↔ y. The initial front-end transformation (2.7) One can check from the explicit formulas (A23) for (α, β, γ) that this transformation maps We have used this plus the fact that (α, β, γ) are symmetric under x ↔ y in deriving (D12d). The other special feature of the front-end transformation (2.7) is that it introduces an overall factor of (1−x) −ǫ . To see what happens to this, focus on the factor [2/ǫ+ln(µ 4 a/E 2 )+ c 1 ] in (D7) for xyȳx. By (D10b), this factor is the same for xȳxy. The front-end transformation (D13) of E together with the overall front-end transformation factor (1−x) −ǫ , followed by the switch of variables x ↔ y, then takes The extra −4 ln(1−y) term above is responsible for the last term in (D11), and we will see later that it conspires in a natural way with similar logarithms in the "(∆t<a) [2]" piece of 2 Re(yxyx) that we will derive from (D8). When the four terms (D12) are added together in (D11), all but the two (iΩ −1,x,1−x ) d/2 terms cancel in pairs. Expanding those in ǫ, we find 2 Re dΓ dx dy ∆t<a [1] virt Ic where Similarly combining (D8), (D10) and fig. 9, we find the remaining contributions from the dimensionally-regulated integration over ∆t < a are 2 Re dΓ dx dy ∆t<a [2] virt Ic yxxy + H with The individual contributions from each diagram to (D17) can be identified by the subscripts.
The phase e −iπ in a logarithm should be interpreted as and the selection of this branch cut is explained in section 4.6 of ref. [23].
We are now in a position to choose D 2 (∆t) of (D1) for the entire sum of crossed virtual diagrams. The 1/ǫ divergence in (D16) represents the dimensional regularization of a ∆t→0 divergent integral a 0 d(∆t)/(∆t). We may use this as a convenient short-cut to read off the ∆t→0 behavior of the integrand from (D16), using the observation of eq. (4.35) of ref [23] that the regulated UV divergence is From the 1/ǫ terms of (D16), we then see that the small ∆t behavior F 2 (∆t) of the d=2 integrand for the sum of virtual crossed diagrams is (Alternatively, one could explicitly extract the ∆t→0 behavior of each diagram and add them up to get the same answer.) Applying the replacement (D5) to (D21) yields our choice of D 2 , given in (A57).
One of the terms we need in our split (D1) of analytic vs. numerical integration is an analytic integral of D 2 (∆t). Integrating (A57) using  From (D2) for 2 Re(xyȳx), combined with (D10) to get other crossed g → ggg diagrams, combined with fig. 9 to relate them to crossed virtual diagrams, we have One can simplify the last −C term. Using (D14) and the fact that every term in the formulas (A13,A14) determining C is proportional to one of (α, β, γ), the last term in (D24) is equivalent to Since we take 2 Re(· · · ) in (D24), we may replace the above by its complex conjugate. In our formalism, conjugating diagrams is equivalent to negating the values of all the x i , and so the conjugate of (D27) is −C(x, z, y, −1, γ, β, α, ∆t). This is the version we have used for our final rewriting of (D24), which is presented as 2 Re B virt Ic in eqs. (A55) and (A56). Following (D1), this (F 2 ) virt Ic = 2 Re B virt Ic is combined with (D 2 ) virt Ic and A pole virt Ic (D23) to give our final total result (A55) for the crossed virtual diagrams.

Sequential Virtual Diagrams
The sum of Class I sequential virtual diagrams (xyxy, xxyy, and xxȳȳ from fig. 4 plus conjugates) are, by fig. 9, just the back-end transformation of the sum of the three g→ggg sequential diagrams shown in the first line of fig. 2 plus conjugates. The latter, computed previously [22], is A seq (x, y) + A seq (x, z) with A seq given by (A59), where the separate terms A seq (x, y) and A seq (x, z) correspond to two different large-N c color routings of the diagrams. 41 The back-end transformation just introduces an overall minus sign, and we must include a loop symmetry factor of 1 2 for the amplitude (blue) or conjugate amplitude (red) loops in the resulting virtual diagrams, giving This result is summarized in (A59). Similarly, as depicted in fig. 10, a front-end transformation of A seq (x, y) + A seq (x, z) followed by x↔y gives the sum 2 Re[ȳxȳx +ȳȳxx + yyxx] of three Class II sequential virtual diagrams: Since rates for all of these processes (as well as for leading-order g → gg) ultimately depend onq and E as q/E, we can rewrite (D29) as 42 This is the result summarized in (A60). A small advantage of (D30) over (D29) for numerical work is that one may work throughout in units whereq=1 and E=1 to get numerical results for rates in units of q/E. Alternatively, one could implement the original (D29) by making E itself an additional argument of all the functions in section A 2 c.
For analytic work, there is a potential conceptual confusion associated with (D30) concerning the meaning of the integration variable ∆t in the definition (A32) of A seq . In all the previous discussion in this paper, ∆t has represented the difference in time between the middle two splitting vertices of interference diagrams like figs. 1, 2, 4, and 5. However, if one steps through in mathematical detail how the explicit formulas of section A 2 c produce equivalence of (D29) and (D30), one finds that ∆t = (1−y) −1/2 ∆t, where ∆t represents the time integration variable associated here with the formula for (D30). In terms of our earlier scaling argument for (D30) from (D29), which bypassed looking at details of the formulas for A seq , this rescaling of the meaning of the ∆t integration variable reflects the fact that formation times scale like E/q.
Finally, we mention that (D14) showed how the combinations (α, β, γ) of helicitydependent DGLAP splitting functions mapped into each other under front-end transformation, but there is no similar relation for the combinations (ᾱ,β,γ) that appear in formulas like (A33) for sequential diagrams. But we have checked that front-end transformation (D13) takes (ᾱ,β,γ) −→ (1−x) 10 (ᾱ,β,γ), where (ᾱ,β,γ) are the combinations of splitting functions that would be obtained by directly evaluating front-end sequential virtual diagrams instead of using our short-cut method of front-end transforming previously known g → ggg sequential diagrams. In detail, (D31) gives (ᾱ,β,γ) in terms of eqs. (A46) for (ᾱ,β,γ) as Unlike (ᾱ,β,γ), the (ᾱ,β,γ) turn out to be symmetric in x↔y and so are unaffected by that step of the transformation of g→ggg diagrams into Class II sequential virtual diagrams in fig. 10.

2 Re(xyyx)
As mentioned in the main text, we can obtain the result for the xyyx diagram of fig. 4 by adapting the results [15] for the similar QED diagram of fig. 12. To go from QED to QCD, we need the following modifications.
• To account for QCD group factors at the vertices, we need to replace However, there are two different large-N c color routings of the QCD diagram, similar to the discussion of color routings of sequential diagrams in section 2.2.1 of ref. [22]. So the overall C 2 A α 2 s corresponds to a factor of 1 2 C 2 A α 2 s per large-N c color routing.
• P e→e and P γ→e are both replaced by P g→gg /C A , where the C A is taken out because each C A in a P g→gg is already explicitly accounted for in the N f α 2 EM → C 2 A α 2 s translation above. Similarly, one should use the (gluonic) QCD formulas of (A46) for the combinations (ᾱ,β,γ) of helicity-dependent splitting functions that are needed for this diagram. (See appendix C for an explanation of why d=2−ǫ versions are not needed.) • Unlike the electron self-energy loop in fig. 12, the corresponding gluon self-energy loop comes with a loop symmetry factor of 1 2 .
• Use the (gluonic) QCD formulas of appendix A 2 for complex frequencies Ω and matrices a of normal modes.
Let A new (x, y) represent a single color routing not including the loop symmetry factor 1 2 . By the same arguments given in section 2.2.1 of ref. [22], the two color routings are related by y ↔ z, and so 2 Re which is (A61). For A new , we can then copy various formulas from ref. [15] with N f α EM → All that remains is the pole piece, which we will package a little differently in this paper than in ref. [15]. Similar to our analysis of crossed virtual diagrams, the pole piece corresponds to the a 0 d(∆t) F d (∆t) + ∞ a D 2 (∆t) term in (D1). The QCD formula we need can be obtained by starting from eq. (F.42) of ref. [15] for QED: 43 where y e ≡ y e /(1−x e ). Here we also adopt the shorthand notation of ref. [15] that In the QED case, it was possible to explicitly perform the integral over y e above to get the pole piece of [dΓ/dx] xyyx . 44 For the QCD translation, the corresponding y integral will be IR divergent because, unlike P γ→e (y), gluon splitting P g→g (y) diverges in the soft limit y→0.
We could do the integral explicitly using our IR regulator δ, but, for our various discussions of cancellations of QCD power-law IR divergences in this paper, it has been very useful to work in terms of dΓ/dx dy for virtual diagrams instead of directly with the IR-regulated dΓ/dx for each diagram. So we'll instead translate the unintegrated version of (D34), from QED to for QCD. Ω i for the QCD diagram is equal to Ω 0 ≡ Ω −1,x,1−x . Now (i) use the fact that the g→gg splitting function P (x) is independent of dimension (see appendix C), (ii) fully expand the above in ǫ and drop terms that vanish as ǫ → 0, and (iii) use y ≡ y/(1−x). This gives Taking 2 Re(· · · ) to add in the conjugate diagram gives what we call A pole new in (A66). 43 Here we have accounted for an overall sign error that appeared in the original published version of eq.

2 Re(xȳȳx)
We obtain 2 Re(xȳȳx) from 2 Re(xyyx) by combined front-end and back-end transformation as depicted in fig. 11. The non-pole (subtracted) piece For the same reasons as described for A seq in (D30), the above can be rewritten as This result appears as the integral in (A69). Here we have called the integration variable ∆t instead of ∆t for reasons similar to those described in section D 3, but here the relation is For the pole piece, we could do the same thing, but we prefer to do the transformation by hand in order to be careful about issues concerning branch cuts. Using fig. 11, the pole piece (D38) for xyyx transforms to The arguments of the above logarithms have minus signs, and we need to decide which branch of the logarithms they land us on. The QED discussion given in appendix H.1 of ref. [15] applies equally well here. That discussion tracks the origin of the complex phases in direct calculations of what we would call here xȳȳx relative to xyyx. The result is that the xȳȳx diagram should have a phase of i d relative to the xyyx diagram, which means ±i −ǫ since the discussion did not keep track of overall signs. Since this means that the branch-cut ambiguity in (D42) resolves as Finally, taking 2 Re(· · · ) of (D44) gives what we callĀ pole new in (A70). There are alternative ways one could write our result for 2 Re(xȳȳx) that may be useful for some purposes. If one wants a formula in terms of the actual duration ∆t of the self-energy bubble, one can make the change of variables (D41) in (A69) to writē Alternatively, going back to (D39), one could use (D31) and scaling arguments similar to (D24-D26) to write (A69) as We have checked that this is the form one would get by directly evaluating the xȳȳx diagram using our methods [15,[21][22][23] instead of taking our shortcut of front-and back-end transforming the xyyx diagram. 2 Re(yxxȳ) +2; 0; +(G− π 2 ) +2; 0; +(G− π 2 ) (A3) (B3) 0; +1; +D 2 Re(yxxȳ) +2; 0; +(G+ π 2 ) ∼ z −5/2 (C1) (C2) 0; +1; +(D+ π 2 ) 2 Re(xȳxy) sum 2 Re(yxxȳ + xȳxy) −2; +1; −X sequential: one color routing of 2 Re(xyxȳ + xxyȳ + xxȳy) Virtual Diagrams (Class I): crossed: sum 2 Re(yxxy +xyxy) +2; −1; +X back-end sequential:  For X(ξ), we have not yet derived an analytic formula. At the moment, we only know that in leading-log approximation for small ≪ ξ ≪ 1 (or, symmetrically, for small ≪ 1−ξ ≪ 1), it is Some individual entries are more divergent than the (small) −3/2 of (E1), but these more severe divergences cancel between pairs of diagrams, leaving behind a net (small) −3/2 diver-gence. For example, the z→0 limit of 2 Re(yxxȳ) and 2 Re(xȳxy) are marked in the table as each diverging like (small) −5/2 , but we give a separate line in the table showing the net  divergence of their sum.  The table explicitly shows as "± mπ n " all contributions that arise from logs of complex phases, which are commented on in section V.
The annotations (A1), (B3), etc. on some entries are just comments to connect to the soft limits of those diagrams considered in previous leading-log analyses of overlap effects. See section E 3 below for an explanation.
2. Assembling y→0 limit of v(x, y) + 1 2 r(x, y) In the table, we have entries for only three of the crossed g→ggg diagrams (plus their conjugates). The full set of crossed g→ggg diagrams ( fig. 1) consists of these three entries plus all possible permutations of the three daughters (x, y, z). But those other cases can be read off from permutations of what is included in the table. For instance, the y→0 limit of 2 Re(xyȳx), which is not listed in the table, corresponds by permutation symmetry to the x→0 limit of 2 Re(yxxȳ), which is listed. We have chosen yxxȳ, yxxȳ, and xȳxy to be our three representative entries in the table in order to highlight their direct back-end relation to the virtual-diagram table entries for yxxy, yxxy, andxyxy: the corresponding rows of the table are just the negative of each other.
The single table entry for g→ggg sequential diagrams shows the A seq of (A32). As discussed in ref. [22], this corresponds to one of two large-N c color routings for the sum of the three diagrams shown explicitly in the top line of fig. 2 (plus their conjugates). The complete set of sequential diagrams and color routings corresponds [22] to summing A seq over all possible permutations of (x, y, z), as made explicit in (A31).
The total differential rate r(x, y) for g → ggg (1.8b) corresponds to the sum over all six permutations of the table entries discussed above. Because of the relationship between limits of those permutations, the y→0 limit of r(x, y) is then twice the sum of the results listed in all three columns y→0, z→0, and x→0 of the subset of g→ggg results given in the table. Adding the g→ggg table entries together then gives Now turn to the virtual diagrams listed in the table. The Class I virtual crossed diagrams in the table correspond to all of the virtual crossed diagrams (top line of fig. 4 plus conjugates) -there are no permutations to add. The Class I and Class II virtual sequential diagrams are related by back-end and front-end transformation to the g→ggg sequential diagrams discussed above. See section D 3 for a discussion. Again there are no permutations to add, and the same is true for the remaining virtual diagram entries 2 Re(xyyx) and 2 Re(xȳȳx).
Because of the addition of "(y ↔ z)" in the definition (1.8a) of v(x, y), the y→0 limit of v(x, y) will sum both the y→0 and z→0 (but not x→0) columns of the virtual diagram which is the negative of (E5). This is in detail how power-law IR divergences cancel in the combination v(x, y) + 1 2 r(x, y) presented in (1.11). Note that we never made use of the x→0 column for the virtual diagrams. Those entries do not add to zero. These divergences (and the related 1−x→0 divergences for class II diagrams) correspond to the blue lines in fig. 13. They do not cause divergences in the applications we have discussed for the reasons described in section III A 4.

The diagrams responsible for double logs in earlier papers
The diagrams that were analyzed in earlier papers [7][8][9] that found the IR double logarithm correspond to the subset of 9 diagrams (A1, A2, ..., C3) 47 depicted by fig. 21, where y represents the softest gluon in the process. Here we comment on why our IR power-law divergences were absent in their analysis.
The y→0 limit of each of these diagrams corresponds to the entries of table I correspondingly marked (A1), (A2), etc. Some entries in the table correspond to more than one of these diagrams: for example, the x→0 limit of 2 Re(yxxȳ) is listed as both (A3) and (B3). That's because permutation symmetries relate this to the y→0 limit of both A3 = xyȳx and B3 = zyȳz. In other places, an entry may be listed as giving only half of the corresponding contribution. For example, the table entries for both the y→0 and z→0 limits of 2 Re(xyyx) are listed as half of the y→0 limit of (A1). That's just a combinatoric issue arising from our 46 The contributions of just Class I diagrams or just Class II diagrams to (E6) are [0; 0; −G+D+X] and [0; 0; + π 4 ] respectively. 47 This naming convention for these diagrams can be made to agree with that used by ref. [7] if our names xE and yE for gluon energies are translated to their zE and ω ′ . In their notation, ref. [7] works mostly in the limit ω ′ ≪ ω ≡ (1−z)E ≪ E.   labeling the two internal lines of the gluon self-energy loop in the xyyx diagram in fig. 4 as y and z = 1−x−y, and in our table there are divergences associated with either becoming soft. In fig. 21, however, y is by definition whichever one of the two is softest. The resulting y→0 divergences for the diagrams of fig. 21 are collected in table II. Each row of table II sums to zero. Consider, for example, the sum A1 + A2 + A3 shown in fig. 22. The reason for this cancellation is that the diagrams are identical except for which line the blue y→0 gluon couples to on the right-hand side, and so the sum is proportional to the sum of those couplings, shown in fig. 23. Because the three hard particles form a color singlet on the right-hand side of this diagram, the coupling of the small-y gluon to the collection of all three will be suppressed compared to its coupling to any individual particle, which is why the leading IR behavior (the power-law divergences) cancel among these diagrams.
In contrast, it's interesting to note that the columns of table II do not sum individually to zero. Consider, for example, the sum A1 + B1 + C1 shown in fig. 24. They differ not only by which line the y→0 gluon couples to on the left-hand side of each diagram but also by whether the y→0 gluon corresponds to a particle propagating in the amplitude (blue line) or conjugate amplitude (red line), which changes the overall time evolution of the diagram. For this reason, one cannot simply factorize out the sum over vertex couplings as we did for A1 + A2 + A3, and so there is no reason for this particular sum of diagrams to be suppressed.
Regardless, the cancellation of each row of table II is sufficient to guarantee that there will be no power-law IR divergences in the sum of all nine diagrams of fig. 21, which is why earlier leading-log analyses did not need to address such divergences.

Derivation of D(ξ)
Here we will give an example of the derivation of one of the boldfaced D's in table I. We will focus on the entry for the x→0 limit of 2 Re(yxxȳ). This is the same, by permutation, as the y→0 limit of 2 Re(xyȳx), to which we now turn since xyȳx is the canonical crossed diagram presented in earlier work [21,23]. Let's look first at the ∆t integral associated with the xyȳx diagram, which is the term ∞ 0 d(∆t) 2 Re C(−1, y, z, x, α, β, γ, ∆t) (E7) of (A11) and (A12), where C is given by (A13) in terms of the D of (A14). One finds that the integral is dominated by ∆t ∼ y for small y. 48 An analytic analysis of the integrand for ∆t ∼ y → 0 yields 49 D(−1, y, z, x, α, β, γ) ≃ D approx (E8) with D approx = − C A α 2 s P (x) 4π 2 y(∆t) 2 ln 48 A quick, initial way to figure out the scaling of the dominant contribution is to make a numerical log-linear plot of ∆t times the integrand vs. ∆t for two extremely small values of y and see how the most prominent feature of the plot scales with y. Because of large round-off error associated with delicate subtractive cancellations in our formulas for small ∆t, we found this method requires using much higher precision numerics than standard machine precision in order to get good results for the integrand at extremely small value of y and ∆t. 49 In particular, D is dominated for ∆t ∼ y by the 2γZ yȳ I 1 and γȲ yȳ Y yȳ I 2 terms of (A14); these (X, Y, Z) are individually given by the 1/∆t terms shown in eq. (D.2) of ref. [21], but the combination X y Xȳ − X 2 yȳ ≃ − x 2 yM0E where Ω 0 = Ω −1,x,1−x as in (2.31). Subtracting the vacuum (q → 0 and so Ω 0 → 0) gives C(−1, y, z, x, , α, β, γ) ≃ C approx with C approx = − C A α 2 s P (x) 4π 2 y(∆t) 2 ln 1 + 2iΩ 0 (1−x) ∆t xy One can rewrite the above as a total divergence, and so do the integral and then take 2 Re(· · · ) to find the leading y→0 behavior of the ∆t integral for 2 Re(xyȳx). This is a y −2 divergence, which would dominate over the y −3/2 divergences of table I except that (E12) exactly cancels the y→0 limit of the pole term for 2 Re(xyȳx). This pole term [23] represents the portion of A pole (A15) attributable to that diagram. The piece of the pole term responsible for the y −2 divergence is 2 Re(· · · ) of the −2γ term in (D8). So we need not worry about the canceling y −2 divergences except that (E10) hides a sub-leading y −3/2 divergence of the integral. (Such cancellations make us wonder whether there is some more elegant analysis of diagrams that would give simpler formulas that more directly reveal the physics of the y→0 limit.) b. The surviving y −3/2 divergence The contributions 2 Re(C − C approx ) to the ∆t integrand that are not accounted for by 2 Re C approx above are dominated 50 by ∆t ∼ y 1/2 . Physically, this corresponds to ∆t ∼ t form (y), where t form (y) is the formation time associated with bremsstrahlung of a soft y gluon.
Repeating the analysis of the small-y expansion of D but now for ∆t ∼ y 1/2 instead of ∆t ∼ y, we find 51 [Ω y csc(Ω y ∆t)] 2 ln(2iΩ 0 ∆t), (E13) 50 One may use the same method as footnote 48. 51 Not much changes from the previous derivation for ∆t ∼ y except that (i) some of the terms that were important for ∆t ∼ y can be ignored for ∆t ∼ y 1/2 , and (ii) it is no longer possible to take the small-∆t approximation to Ω + csc(Ω + ∆t) when calculating Z yȳ . In particular, we find that Ω + is of order the inverse y-formation time for small y, so that Ω + ∆t ≪ 1 for the previous case ∆t ∼ y but Ω + ∆t ∼ 1 for the ∆t ∼ y 1/2 case here. This point only matters for Z yȳ since we find that the small-y limits of the relevant (X, Y )'s are not sensitive to Ω + csc(Ω + ∆t).
where Ω y ≡ −iq A 2yE . (E14) Comparing to the already-accounted-for D approx of (E9), and remembering that now ∆t ∼ y 1/2 , D ≃ D approx + δC (E15) with δC ≡ − C A α 2 s P (x) 4π 2 y Ω y csc(Ω y ∆t) 2 − 1 (∆t) 2 ln(2iΩ 0 ∆t). (E16) We've called it δC instead of δD because it already vanishes in the vacuum limitq→0, which takes both Ω y and Ω 0 above to zero. So the vacuum subtraction has no effect on this contribution to D. The y −3/2 divergence of 2 Re(xyȳx) will now come from taking the integral over ∆t of 2 Re δC. By changing integration variable to τ ≡ iΩ y ∆t, which runs from 0 to e iπ/4 ∞, and then arguing that one can safely add a contour at infinity to deform the integral to be from 0 to +∞, one gets The integral formula 52 Re iΩ y ln 2πΩ 0 Ω y − γ E .
In the style of (E1), this is with D(ξ) determined in this derivation to be (E2). Permuting x ↔ y in (E20) gives the entry in table I for 2 Re(yxxȳ) as x → 0.

Derivation of G(ξ)
Now we give an example of the derivation of one of the boldfaced G's in table I. We focus on the entry for the y→0 limit of 2 Re(yxxȳ), which by permutation is the x→0 limit of the same canonical crossed diagram 2 Re(xyȳx) analyzed in the previous subsection. a. Spurious x −5/2 divergence of 2 Re(xyȳx) Similar to the y→0 limit of 2 Re(xyȳx) studied in section E 4, the ∆t integral (E7) also generates a spurious dominant divergence in the x→0 limit. In this case, the integral is dominated by ∆t ∼ x 3/2 , for which 53 D approx = − C A α 2 s P (y) 4π 2 x(∆t) 2 ln xy 1−y + 2iΩ x ∆t + 1 + 2iΩ x (1−y) ∆t xy where is the small-x limit of Ω 0 . Correspondingly, When integrated, this generates an x −5/2 contribution to the ∆t integral, which is canceled by a similar contribution from the pole term. The relevant piece of the pole term again comes from the −2γ term in (D8).
b. The surviving x −3/2 divergence In this case, the dominant contribution to 2 Re[C − C approx ] comes from two places. One is ∆t ∼ x 1/2 , which physically corresponds to ∆t ∼ t form (x). The other is sub-leading corrections to the ∆t ∼ x 3/2 region we just analyzed above.
Let's start with ∆t ∼ x 1/2 . In this region, we find The difference of this with the already-accounted-for D approx of (E21) is Taking ∆t → 0 above, 2 Re(δC) diverges as 2 Re δC (∆t∼x 1/2 ) ≈ C A α 2 s P (y) 2π 2 x ∆t Re(iΩ x ), (E26) 53 The situation is similar to footnote 49 except that here X y Xȳ − X 2 yȳ ≃ − x 3 y(1−y)E 2 (∆t) 2 1 + 2iΩx(1−y) ∆t xy and γ ≃ 2P (y)/y 2 (1−y) 3 x 3 C A . and so we cannot simply integrate 2 Re δC to find the result we are interested in. In general, the 1/∆t divergence of individual diagrams is what created the need for analyzing what we call pole terms of diagrams. In the current case, this divergence shows up at an order in x that makes it relevant to the integral of 2 Re δC. We will need to subtract out the 1/∆t divergence to get a convergent integral and then add the subtraction back in as part of the pole term, as in (D1). Following (D5), at this order in x we will choose whose ∆t → 0 behavior matches (E26). Eq. (E27) is the same as taking the small-x limit of applying (D5) to the more general small-∆t result (D3) for 2 Re(xyȳx). Defining τ ≡ iΩ x ∆t, the integral we want is 2 Re dΓ dx dy (See appendix B for the last integral.) Now turn back to ∆t ∼ x 3/2 . Carrying out the expansion of D to next order in x (including the size of ∆t in the counting of order), we find 54 Note that R 0 is O(1), but δR andξ are O(x) and so small. We could have more thoroughly written out the x expansion of what is shown explicitly in (E29), but keeping it in its current form will be convenient. For example, not explicitly expanding γ (A23) will make it simpler to see what parts of this calculation eventually cancel with the pole terms at this order in 54 We will not list intermediate steps here except to mention, as a checkpoint, that X y Xȳ − X 2 yȳ = − x 3 yzE 3 (∆t) 2 1 + (2 +ξ −1 )iΩ x ∆t −ξ −1 (Ω x ∆t) 2 1 + O(x 2 ) , which at leading order in x matches the simpler formula of footnote 53.
x. Subtracting the vacuum limit from (E29) gives At leading order in x, this reproduces (E23), but (E31) correctly accounts for the next order in x as well. At that order, D 2 (E27) is relevant, and its subtraction must be included as well, in order for the ∆t→0 integration to converge. It's convenient to use the leading-order conversion 55 to rewrite (E27) as [The leading-order conversion is adequate because D 2 is already a sub-leading effect to our calculation of d(∆t) 2 Re(C − D 2 ).] For ∆t ∼ x 3/2 , the argument of the csc is small, so we may approximate This matches the 1/∆t divergent behavior of 2 Re(C (∆x∼x 3/2 ) ), as D 2 should. It is also convenient to switch from the ∆t variable, which is O(x 3/2 ) in (E31), to the O(1) variablē in terms of which Expansion in x is now equivalent to expansion inξ. Expanding explicitly to NLO inξ, we find that we can rewrite the argument of Re above as iΩ x ξ dτ × d dτ − (1 +ξ) τ ln(1 +τ ) +ξτ (1 +τ ) .
Integration then gives 2 Re dΓ dx dy The last element we need is to extend analysis of the O(x −5/2 ) pole terms to O(x −3/2 ). Since we have had to make the D 2 subtraction above, we also need to add the D 2 term back to the pole terms as in (D1). Using (E33) and 2 Re(· · · ) of (D6-D8), and expanding in x, we find lim "a→0" 2 Re dΓ dx dy Since the O(x −5/2 ) pieces have canceled, we may now use leading-order expressions for z and γ to get 2 Re dΓ dx dy xyyx ≃ C A α 2 s P (y) 2π 2 x Re 2 ǫ + 2 ln µ 2 iΩ x E − ln e −iπ xy(1−y) In the style of (E1), this is with G(ξ) determined in this derivation to be (E3). Permuting x ↔ y in (E42) gives the entry in table I for 2 Re(yxxȳ) as y → 0. This has been a complicated derivation of G(ξ). Reassuringly, one can confirm the final answer numerically by comparing to the soft limit of our full numerical results for the diagram. and S NLO = 1 0 dx 1/2 0 dy v(1−x, y) θ(y < x 2 ) + 1 2 r(x, y) θ(y < x) θ(y < 1−x 2 )