In-medium loop corrections and longitudinally polarized gauge bosons in high-energy showers

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. In this particular paper, we show (i) how the “instantaneous” interactions of Light-Cone Perturbation Theory must be included in the calculation to account for effects of longitudinally-polarized gauge bosons in intermediate states, and (ii) how to compute virtual corrections to LPM emission rates, which will be necessary in order to make infrared-safe calculations of the characteristics of in-medium QCD showering of high-energy partons. In order to develop these topics in as simple a context as possible, we will focus in the current paper not on QCD but on large-Nf QED, where Nf is the number of electron flavors.

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]. A longstanding problem in field theory has been to understand how to implement this effect in cases where the formation times of two consecutive splittings overlap.
Let x and y be the longitudinal momentum fractions of two consecutive bremsstrahlung gauge bosons. In the limit y ≪ x ≪ 1, the problem of overlapping formation times has been analyzed at leading logarithm order in refs. [3][4][5] in the context of energy loss of highmomentum partons traversing a QCD medium (such as a quark-gluon plasma). Together with Chang, we subsequently developed and implemented field theory formalism needed for the more general case where x and y are arbitrary [6][7][8][9]. We used the formalism to calculate the LPM interference effects on real double gluon bremsstrahlung g → ggg in medium, given by the processes depicted in fig. 1, in the high-energy limit. [For the sake of simplicity, that specific calculation also made other simplifying approximations by taking the multiple scattering (q) limit and the large-N c limit.] But that calculation was incomplete for two reasons.
First, the calculations of refs. [6][7][8][9] only included gluons that were transversely polarized. But the polarization of the intermediate-state gluon in the first diagram of fig. 1 does not have to be transverse: a full calculation should include the effects of the longitudinal polarization as well. 1 It has long been known how to integrate out the effects of longitudinal polarizations and find the effective theory of the transverse polarizations, which is the program of Light Cone Perturbation Theory (LCPT) [10][11][12]. 2 In LCPT, where light-cone gauge A + = 0 is used, the longitudinal polarization gives rise to an interaction that is instantaneous in light-cone time x + ≡ x 0 +x z (similar to how in Coulomb gauge the non-transverse 1 Some historical explanation: Longitudinal polarizations were left out of our earlier analysis of refs. [6][7][8][9] because of an unstated and incorrect assumption. Consider the two consecutive vertices appearing in the first diagram of fig. 1. Our interest is in the case where the formation times associated with these vertices overlap, in which case the characteristic scale of the time separation of those vertices (and so also their spatial separation in the direction of motion z) is of order the formation time. That formation time is parametrically large in the large-energy limit (t form ∼ E/q, ignoring dependence on daughter momentum fractions). The effects of longitudinal polarizations do not extend over large distances, and so ref. [6] implicitly ignored them. But that's not a valid argument that the effects of nearly-coincident gluon emissions can be ignored, as was made clear in our own calculations by the subsequent calculation in ref. [9] of the 4-point vertex contributions represented by the last diagram of fig. 1. It was not until we started computing virtual corrections (discussed in a moment) that we clearly realized our mistake in ignoring intermediate longitudinal polarizations. 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 [13]. For some less-simple examples, see also ref. [14].  [6][7][8][9]. Only the high energy particles are shown: their many interactions with the medium are not shown explicitly. x and y represent the momentum fractions (relative to the initial particle) of two of the daughters of this splitting process.
FIG. 2: The contribution to g → ggg due to longitudinal polarization of the intermediate gluon line (represented by a bar across the line). This process must be added to the g → ggg amplitude of fig. 1 if the intermediate lines in fig. 1 are summed only over transverse polarizations as in refs. [6][7][8][9].
polarization A 0 of the gauge field responds to sources instantaneously in normal time x 0 , even though physical quantities ultimately do not). In LCPT, the longitudinal polarization can be integrated out to yield a Hamiltonian that is supplemented by just an extra 4-point interaction among the fields. The result is to add one extra diagram (and its permutations), fig. 2, to the sum of diagrams depicting the amplitude in fig. 1. In fig. 2, the bar across the intermediate gluon line indicates that it is an (integrated-out) longitudinally polarized gluon. That gluon line is drawn vertically to emphasize that it represents an instantaneous interaction in the context of time(x + )-ordered perturbation theory, with x + running from left to right in the figure. One goal of this paper is to include this type of process into our earlier calculations of the effect of overlapping formation times on sequential bremsstrahlung. A calculation of in-medium real double bremsstrahlung g → ggg is not enough, however, to study the effects on in-medium energy loss and the characteristics of in-medium shower development. Crudely analogous to what happens for vacuum bremsstrahlung in QED, there are infrared divergences whose treatment at this order requires additionally computing virtual corrections to single bremsstrahlung g → gg, such as depicted by the diagrams in fig. 3. The second goal of this paper is to show how to calculate such loop corrections (without soft gluon approximations) in the context of finding the effects of overlapping formation times on energy loss and shower development. A new technical feature, compared to our calculations of g → ggg, is that loop calculations require UV regularization and renormalization. In particular, one of the effects of loops (besides mitigating certain infrared divergences of the g→ggg calculation) will be to replace the coupling α s in the leading- order single-bremsstrahlung process g → gg by the running coupling α s (Q ⊥ ) evaluated at a characteristic transverse momentum scale of the LPM-modified bremsstrahlung process.
In this particular paper, we work out all of these issues for a warm-up theory: QED in the large-N f limit, where N f is the number of electron flavors. The large-N f limit is taken just to reduce the complexity of the calculations and so streamline the presentation of the important developments. The method itself should be straightforwardly generalizable to the QCD analysis of the LPM effect in figs. 1-3 (at least in the large-N c limit of refs. [6][7][8][9]), but we leave that QCD calculation for the future.
There are several aspects of this paper that are related to recent work by Beuf [19,20] and Hänninen, Lappi, and Paatelainen [21,22] on next-to-leading-order deep inelastic scattering (NLO DIS). Here, making use of LCPT, we will study the combination of 1→3 splitting e → eēe and UV-renormalized one-loop corrections to 1→2 splitting (e → γe or γ → eē) in the presence of a thick medium (though the formalism we use can in principle handle more general situations). In NLO DIS, Hänninen et al. instead used LCPT to study the combination of 1→3 splitting γ * → qqg and UV-renormalized one-loop corrections to 1→2 splitting (γ * → qq) in the presence of an extremely thin medium. We briefly comment further on the similarities and dissimilarities in appendix A.
The goal of this paper is to develop calculational methods. We leave discussion of application of the results to computing IR-safe characteristics of energy loss and shower development to a later paper [15].
B. Overlapping formation times in large-N f QED Fig. 4 depicts the showering of a high-energy particle in a medium. For simplicity of discussion, let's focus for the moment on roughly-democratic splitting processes, meaning that neither daughter of any splitting is soft (i.e. neither has very small energy compared to the parent). Let's also first discuss the case of QED with only one electron flavor (N f = 1). If the coupling associated with the splitting vertices is small, then there is a hierarchy of scale in fig. 4 between (i) the formation length l form associated with each splitting (denoted by the length of the ovals) and (ii) the mean free path l rad between consecutive roughly-democratic splittings. Parametrically, l rad ∼ l form α . (1.1) A rough mnemonic for this result is that each formation length of media traversed offers one opportunity for splitting, with probability of order α. On average it then takes ∼ 1/α such opportunities to radiate. Because of this scale hierarchy, the probability for two consecutive splittings to have In (a), splittings are independent and do not interfere with each other. In (b), two splittings have overlapping formation lengths and so may not be treated as independent. Showering is extremely collinear in the high energy limit, but the transverse directions in these schematic drawings have been drastically magnified for the purpose of illustration.
overlapping formation times, as depicted in fig. 4b, is suppressed by a factor of α. To leading order in α, one may treat the splittings in showers as non-overlapping, as depicted in fig. 4a, and so may treat the probabilities of each splitting in the shower as independent. The goal of the current program is to understand how to calculate (beyond soft-emission approximations) the α-suppressed effect of overlapping formation times, such as the overlapping pair of splittings in fig. 4b. So far, the qualitative discussion has been the same as the discussion for QCD in the introduction of ref. [6]. But now consider the large-N f limit, which will reduce the number of calculations needed in this paper. If the properties of the medium are held constant, then N f does not affect the rate for bremsstrahlung e → γe. Formation times are also not affected. But the rate for pair production γ → eē is proportional to N f , and so the mean free path for pair production is smaller than for bremsstrahlung by a factor of N f . The hierarchy of scales relevant to typical showering is then summarized by fig. 5, assuming that N f is large but N f α is still small: 3 α ≪ N f α ≪ 1. (1. 2) The probability that two of the closer splittings in the figure (e → γe → eēe) might overlap will be of order N f α. That's overlap of (i) bremsstrahlung with (ii) the subsequent pair production later initiated by the bremsstrahlung photon. The chance of any other type of overlap, such as two consecutive bremsstrahlung processes e → γe → γγe, is only of order α and so is parametrically less likely in the large-N f limit. So, to get the dominant correction due to overlapping formation lengths in this problem under the formal assumption (1.2), the only type of overlap we need to compute is the type shown in fig. 6 (plus related virtual processes). The interference diagrams needed to calculate the effect of overlapping formation times in large-N f QED are shown in figs. 7-9. These diagrams are drawn using the conventions of  refs. [6,7] and represent contributions to the rate for the process. Each diagram shows the product of a term in the amplitude for the process (colored blue in the diagram) times a term in the conjugate amplitude (colored red). The diagrams are time-ordered (more precisely in this paper, light-cone time ordered), with time running from left to right. We will only consider rates which are integrated over the transverse momenta of the final particles, which allows one to ignore the evolution of any final-state particle after it has been emitted in both the amplitude and conjugate amplitude. 4 So, for instance, the diagrams of fig. 7a and 7g represent the (p ⊥ -integrated) interferences of fig. 10.
We will see later that the results for most of the virtual diagrams can be related to results for e → eēe ( fig. 7), or to each other. There will only be one quintessential virtual diagram that we will have to compute from scratch, which will turn out to be the boxed diagram shown in fig. 8k. That diagram, and the related diagram of fig. 9r, contain the only true Diagrams with only transverse intermediate photons: Diagrams including longitudinal intermediate photons: FIG. 7: Time-ordered interference diagrams for e → eēe in large-N f QED. As in refs. [6,7], blue represents a contribution to the amplitude and red represents a contribution to the conjugate amplitude. As in the other figures of this paper, repeated interactions with the medium are present but not explicitly shown. The complex conjugates of the above interference diagrams should also be included by taking 2 Re[· · · ] of the above.
FIG. 8: Time-ordered interference diagrams for the virtual correction to e → γe in large-N f QED. The boxed diagram is the only ones whose result cannot be simply related to one of the e → eēe diagrams of fig. 7. Again, complex conjugates should be included by taking 2 Re[· · · ].
UV divergence in these large-N f calculations 5 and will be the diagrams responsible for the usual renormalization of the QED coupling α.
In the approach of our earlier work [6][7][8][9], inspired by Zakharov's treatment of the LPM effect [18], we interpret the interference diagram of fig. 7a, for instance, as the evolution of an initial eē pair, where the e represents the initial electron in the amplitude and theē represents the same particle in the conjugate amplitude. In this language, the time evolution of fig. 7a is interpreted, as shown in fig. 11a, as a phase of 3-particle evolution (γe in the amplitude andē in the conjugate amplitude), followed by 4-particle evolution (eēe in the amplitude andē in the conjugate amplitude), followed by 3-particle evolution (eē in the amplitude and FIG. 9: Time-ordered interference diagrams for the virtual correction to γ → eē in large-N f QED. All of these diagrams can be related to one of the diagrams of figs. 7 or 8. γ in the conjugate amplitude). We then used symmetries of the problem to replace each medium-averaged N-particle evolution problem by an effectively (N−2)-particle evolution problem, 6 as summarized in fig. 11b. One of our tasks in the current paper will be to translate the LCPT diagrammatic rules for longitudinal photon interactions into corresponding rules within our framework.
Before going forward, we should mention that formally there are some additional virtual diagrams in Light Cone Perturbation Theory, which will be shown and discussed later (fig. 23) but which are negligible in the high-energy limit provided we use dimensional regularization.

C. Assumptions and Approximations
As with the earlier work in refs. [6][7][8][9], the formalism we will present is quite general, but for now we will implement it in simple cases that make the actual calculations much easier.
First, we will assume that the medium is large, static, and uniform. More precisely, we approximate the medium as (statistically) uniform and unchanging over the scale of formation lengths/times. 1-particle 1-particle 2-particle effectively effectively effectively 3-particle 3-particle Second, we will make the "multiple scattering approximation" that interactions with the medium can be characterized by the parameterq that is often used in discussions of quarkgluon plasmas, and which is the proportionality constant in the relationship Q 2 ⊥ =q ∆z, where Q ⊥ is the transverse momentum that a high-energy charged particle picks up passing through length ∆z of the medium. For a variety of reasons, the effective value ofq can have logarithmic dependence on energy at fixed order in the coupling α that controls high-energy splittings. There are important qualitative differences between QCD and QED, which we briefly review in appendix C, but it is off the main topic of this paper. Here, we will make the approximation thatq is a constant and assume that one is using a value ofq appropriate for the overall energy scale of the initial particle.

D. Some Qualitative Results
This paper further develops methods for next-to-leading order calculations of the LPM effect and gives analytic results for the case of large-N f QED. The details of the full analytic results are complicated enough that we will not present them here in the introduction. And we are mostly leaving numerical analysis and application of those results to a later paper. But there are two features of our results which we will present here.
The first feature regards the appropriate choice of renormalization scale for the factor of α that controls the cost of high-energy splitting and determines the overall importance of overlapping formation times in showering-that is, the scale of the α associated with each high-energy splitting vertex. In earlier work [6], we asserted based on qualitative physical reasoning that this coupling should be taken to be α(Q ⊥ ), where Q ⊥ is the characteristic transverse momentum scale of the LPM-modified splitting process. Our explicit next-toleading-order LPM results bear this out. For roughly-democratic splitting processes, we find that the logarithmic terms in our NLO results for single splitting are 7 where µ is the renormalization scale for α(µ), β 0 is the first coefficient of the renormalization group β-function for α, and |Ω| is a frequency that is of order the inverse formation time. To avoid poor convergence of the perturbative expansion due to large logarithms, one should therefore choose µ ∼ |EΩ| 1/2 above. As we review in section II, this is indeed equivalent (for roughly-democratic splittings) to µ ∼ Q ⊥ . The other feature regards our result for the effect of overlapping formation times on the real double splitting process e → γe → eēe. Though our analytic result is complicated for the general case, there is a simple leading-log formula in the limit that the intermediate photon is soft. That formula can be derived with a relatively simple analysis based on (i) the leadingorder LPM formula for pair production γ → eē combined with (ii) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution of parton distributions to get the probability of seeing the photon γ inside the initial electron e. We explain this analysis in section II, which yields the formula where x γ is the longitudinal momentum fraction of the intermediate photon relative to the initial electron in the bremsstrahlung e → γe, and y e is the longitudinal momentum fraction of the final electron relative to the photon in the subsequent pair production γ → eē. The ≈ sign in (1.4) is our notation for indicating that this formula is only valid at leading-log order. The ∆ in front of dΓ/dx γ dy e is our notation [7] to denote the effect of overlapping formation times on rates for double splitting. That is, ∆[dΓ/dx γ dy e ] is the difference between (i) dΓ/dx γ dy e and (ii) what one would get by always treating the two consecutive mediuminduced splittings e → γe and γ → eē as quantum-mechanically independent. Fig. 12 verifies the leading-log approximation by showing that the ratio of our full numerical results for e → eēe divided by the approximation (1.4) extrapolates to 1 as x γ →0. The convergence is slow because logarithms grow slowly. As we will discuss in section II, there are two interesting qualitative differences between QED e → eēe and QCD g → ggg: (i) the x γ ≪1 behavior (1.4) will not give rise to any infrared (IR) divergences in energy-loss calculations, and so needs no IR cancellation between real and virtual diagrams, and (ii) the in-medium collinear logarithmic enhancement factor in (1.4) does not cancel between real double-splitting diagrams as it did for g → ggg [7] via a Gunion-Bertsch cancellation. 8 This reflects a difference between soft pair production and soft bremsstrahlung. 7 Specifically, see eq. (4.38) later in this paper, specializing here to the roughly-democratic case where neither x e nor 1−x e are small. We have not drawn from our results any conclusions about the best renormalization scale in limiting cases such as 1−x e ≪ 1, because then the logarithmic corrections due to choice of renormalization scale are overwhelmed by other, power-law corrections which dominate. Those power-law corrections are the same parametric size as (2.21). 8 Specifically, see appendix B of ref. [7]. The vertical axis is the ratio of full numerical results to the leading-log approximation, and the horizontal axis is the inverse logarithm 1/ ln(x −1 γ ). The data points come from numerical calculations of our general e → eēe results, given in appendix E. The extrapolation to x γ →0 is taken by a straight-line fit to the leftmost two data points of each set. Extrapolations are shown for two different values of y e . Even though QED is qualitatively different from QCD in these respects, it nonetheless provides a good training ground for working out calculational methods for overlap effects in QCD. And overlap effects in QED are interesting in their own right, even if so far only treated here in the large-N f limit.

E. Outline
We have tried to organize this paper so that the main text gives an introduction and overview of the techniques we use while the fine details of the calculation are left to an extensive set of appendices. In section II, we summarize qualitative differences between QCD and (large-N f ) QED regarding the LPM effect in both single-splitting and overlapping double-splitting processes. A quantitative review of single-splitting formulas is left to appendix B, and appendix C contains a translation of modern notation usingq to the QED results originally presented by Migdal [2]. Section III introduces the elements needed for organizing calculations of overlapping formation-time effects in terms of light-cone perturbation theory (LCPT), with the full list of details left to appendix D. The use of those rules to calculate overlap effects in the real double-splitting process e→eēe of fig. 7 is left to appendix E, where much of the calculation is adapted from our previous work on g→ggg in QCD. In section IV, we turn to the virtual corrections of figs. 8 and 9. We first discuss the techniques needed, which we call back-and front-end transformations, to easily relate almost all virtual diagrams to non-virtual diagrams. We then turn to the one remaining virtual diagram (the boxed diagram of fig. 8k) and outline its computation and renormalization but leave details for appendices F and G. 9 Finally, we end with a brief conclusion in section V.

F. Reference Acronyms
When discussing detailed formulas in appendices and sometime footnotes, we will often need to refer to particular sections or equations of our earlier work [6][7][8][9]. To streamline such references, we will often refer to our earlier work in such cases by the author acronyms and numbers AI1 [6], ACI2 [7], ACI3 [8], and ACI4 [9]. So, for instance, "ACI3 (3.4)" will be shorthand for eq. (3.4) of ref. [8], and "AI1 section II.A" will be shorthand for section II.A of ref. [6].

A. Single splitting
The crucial difference to keep in mind between the LPM effect in QCD and QED is that the photon is neutral but the gluon has color charge. In terms of parametric behavior, this makes a huge difference in how LPM suppression behaves when a bremsstrahlung photon or gluon is soft. The LPM effect relies on the near-collinearity of high-energy splitting processes, and the formation time l form (and so the amount of LPM rate suppression) is smaller when the splitting is less collinear. Since it is easier to deflect a low-momentum particle than a high-momentum particle, the collinearity of QCD gluon bremsstrahlung is controlled by how soft the bremsstrahlung gluon is: high-energy gluon bremsstrahlung is less LPM-suppressed the softer the bremsstrahlung gluon is. A soft photon, however, does not scatter in first approximation, and so the collinearity of QED photon bremsstrahlung is insensitive to this. But a softer photon means a longer-wavelength photon, which means a photon with less resolving power, which is why the QED case has the opposite behavior: high-energy photon bremsstrahlung is more LPM-suppressed the softer the photon is.

Formation Times
A little more concretely, let's review, in a general way, some standard parametric estimates of formation times. Consider a single splitting process of a high-energy particle with energy E into daughters of energy E 2 =x 2 E and E 3 =x 3 E and ask for the formation time. That's the time scale |t −t| over which it's possible for splitting in the amplitude at time t to interfere with splitting in the conjugate amplitude at timet, as depicted by the interference diagram in fig. 13. That time scale t form is determined by the difference δE of the energies of (i) the two-particle state |2, 3 after splitting in the amplitude and (ii) the one-particle state before splitting in the conjugate amplitude: where p i are the transverse momenta, which we take to be small compared the energies E i . In this paper, we will assume throughout that energies are high enough that we can ignore the effective masses m i of the high-energy particles in the medium. One can then make a rough parametric estimate t ∼ 1/|δE| of the formation time by using the definition ofq to take p 2 i ∼q i t above, whereq i is theq appropriate for particles of type i. Then, from (2.1), and so .
For photon bremsstrahlung e → γe, we haveq γ =0, and (2.3) gives whereq ≡q e . Above, x e and x γ (which in this case is x γ = 1−x e ) are respectively the longitudinal momentum fractions of the electron and photon daughters of the splitting. In contrast, for gluon bremsstrahlung q → qg or g → gg, we haveq q ∼q g and The parametric estimates (2.4) and (2.5) show that the photon formation time becomes long (more LPM suppression) as x γ → 0 whereas the gluon formation time becomes short (less LPM suppression) as x g → 0.
We've used the very general formula (2.3) to derive standard parametric estimates of LPM formation times because it is then easy in the same breath to obtain the standard parametric result for QED pair production γ → eē, which is qualitatively different from photon bremsstrahlung. (2.3) gives The QCD pair production formation time (g → qq) is parametrically similar to this, as well as to the the QCD bremsstrahlung formation time (2.5).

Splitting Rates
Below, we will need the differential rates dΓ/dx for single splitting. For nearly-democratic splittings, the relation given by (1.1) is that Γ ∼ α/t form . But we will also be interested in non-democratic splittings (one daughter soft) and will want to discuss the differential rate dΓ/dx, where x is one of the daughter momentum fractions. The differential rate is sensitive to the DGLAP splitting function for the process: for bremsstrahlung (e → γe or X → gX, with x g being the softest gluon in the case g → gg) and dΓ dx pair ∼ N f α t form (2.9) for pair production (γ → eē or g → qq). For a review of precise, quantitative formulas for single splitting rates (in the high energy limit where theq approximation can be used), see appendix B. Appendix C discusses exactly how, in the QED case, these formulas match up to the original results presented by Migdal [2] in the case where the medium is an atomic gas.

An aside on nearly-democratic splitting
We mention in passing that QED and QCD are parametrically similar in the case of nearly-democratic splitting (neither daughter soft compared to the parent). All the formation times discussed above then have the same order of magnitude: (2.10) (The splitting rates are then also all parametrically the same except for factors of N f .) The amount of transverse momentum Q ⊥ transferred from the medium in one formation time is of order Q 2 ⊥ ∼qt form , which can be written in a number of equivalent ways: (nearly-democratic splitting). (2.11) We used the last form to identify µ ∼ |EΩ| 1/2 with Q ⊥ in the discussion of renormalization scale in section I D.
One could of course also write down case-by-case parametric estimates of Q ⊥ for nondemocratic splittings, but we do not have need of them.

B. Overlapping double splitting
In ref. [7] (ACI2), 10 we discussed how to parametrically estimate the size (though not the sign) of QCD overlapping formation-time effects for real double splitting g → ggg by using formation times to estimate rates. Here we will first review that QCD estimate and then see what changes when we switch to (large-N f ) QED.

Review of QCD estimate for g → ggg
Let yE and xE be the energies of the softest and next-softest of the three final-state gluons, so that y x 1−x−y. Since this is QCD, the formation time t form,y for radiating the y gluon will be shorter than the formation time t form,x for radiating the x gluon. The probability that a y emission happens to take place during the x emission is then just t form,x times the rate of y emission. Overall, that means that the joint differential rate for overlapping x and y emissions is (2.14) The fact that two emissions overlap does not a priori mean they will influence each other, and so the qualitative argument for (2.14) only provides an estimate for how large overlap effects might be. However, detailed calculations [7] of overlap effects for g → ggg confirm it. Note that if one multiples (2.14) by the energy (x + y)E lost by the leading parton and integrates over x and y, one would find a power-law divergent contribution from overlap effects to energy loss. No such power-law divergence appears in the soft-gluon bremsstrahlung calculations of ref. [3][4][5], but those calculations inextricably combine the effects of soft virtual emission with those of soft real emission, for which there are cancellations in QCD.
Here we have only estimated the size of soft real emission alone.

QED estimate for e → eēe
For large-N f QED we've already identified that the dominant overlap correction comes from the e → eēe process of fig. 6. Let x e E and y e E be the energies of the two final-state electrons. As shown in the figure, our convention throughout this paper will be that x e is the daughter whose electron line is connected to the initial-state electron, and y e is the electron in the eē pair produced by the intermediate photon. This is a distinction made possible by the large-N f limit, in which the chance that the electron in the γ → eē pair has the same flavor as the initial electron (which would allow interference terms involving exchange of the two electrons) is 1/N f suppressed. In this section, we will focus on the case x e ∼ 1, which includes the case x γ → 0 but not x e → 0. (Our later explicit calculations of diagrams make no such assumption.) The formation time for the initial bremsstrahlung e → γe can be taken from (2.4): To get the formation time for the subsequent pair production process γ → eē, we have to be careful applying the single splitting estimate (2.6) because we are considering the case where the photon is already soft. The E in (2.6) is the photon energy, which we are now calling x γ E. The x e in (2.6) is the longitudinal momentum fraction of the pair's electron relative to the photon, which in the application here (see fig. 6) is With these substitutions, (2.6) gives t form,y ∼ y e (1−y e )x γ Ê q (2.17) in the present context. As in the QCD estimate, this is t form,x , and so we can estimate the probability of overlap the same way, using (2.12), which can also be written Here we need to make sure to use the pair-production formula (2.9) to get dΓ/dy e ∼ N f α/t form,y . So, using (2.17) above, y e x e FIG. 14: The double splitting process e → γe → eēe viewed in an approximation (relevant to overlap effects in the limit x γ ≪ 1) treating it as fundamentally pair production γ → eē with the γ interpreted as a parton inside the original electron, described by DGLAP evolution.
Unlike the QCD case, the integral over y e does not diverge: Furthermore, were we to use this formula to calculate the overlap effects on energy loss, we would get finite results, even without accounting (as one still should) for virtual corrections. QED has much better infrared behavior than QCD for these calculation because the LPM suppression of photon radiation increases as the photon becomes softer.

Leading-log formula for e → eēe
Before moving on to the details of our complete calculation of overlap effects in the general case, we first briefly discuss a relatively simple argument for the leading-log formula (1.4) for x γ ≪ 1. As mentioned earlier, the existence of such a logarithm is another difference between e → eēe and g → ggg.
Consider the e → eēe process of fig. 6 and re-interpret it as the Feynman diagram shown in fig. 14. In the latter figure, we have not explicitly drawn the (still overlapping) ovals corresponding to formation times. Instead, the box here denotes the leading-order pairproduction process γ → eē via interaction with the medium, and we will consider that photon γ as a parton contained inside the initial electron. In a related context, ref. [7] (ACI2) argues that the corresponding parton-model-like (i.e. DGLAP) approximation to such a process has the generic form 11 where dΓ/dy is the differential rate for the fundamental parton-level process (represented here by the box in fig. 14) and here the rest of (2.22) represents the probability of finding the photon inside the initial electron. The analysis of that probability is different from what it would be in vacuum because scattering from the medium cuts off collinear logarithms.
[See ref. [7] for the argument that the DGLAP logarithm ln(Q 2 /Q 2 0 ) of the virtuality ratio translates to the logarithm ln(t form,x /t form,y ) of the formation-time ratio.] The translation of (2.22) to fig. 14  The parametric estimate of dΓ/dy e was given earlier as ∼ N f α/t form,y ∼ N f α q/y e (1−y e )E, but here we want the exact leading-order result. That's reviewed in appendix B and (using the same translations as in our parametric analysis) is Our earlier work did not include the effects of longitudinally-polarized gauge bosons in intermediate states. One may retain a description solely in terms of transverse polarizations by first integrating out the longitudinal polarizations in light-cone gauge, as in Light Cone Perturbation Theory. This gives rise to 4-fermion interactions such as shown in fig. 15a, which are instantaneous in light-cone time x + , local in transverse position x ⊥ , and non-local in x − . In Hamiltonian formalism (as in LCPT), the legs can be incoming or outgoing, such as shown in fig. 16 for fig. 15a. In LCPT, a similar 4-particle interaction arises from fermion exchange, corresponding to fig. 15b. Loosely speaking, this interaction term accounts for the difference between using (i) the actual off-shell exchanged fermion and (ii) treating that fermion as though it had on-shell polarization u p or v p . An important feature of LCPT is that fig. 15 for QED (and similar 4-field diagrams for QCD) is all that is needed to account for the effects of longitudinal gauge bosons or off-shell fermion polarizations: there are no n-field interaction terms in the light-cone Hamiltonian with n > 4.
Figs. 17 and 18 show two examples of standard LCPT rules for vertices (shown inside the boxes) and examples of their translation to corresponding vertices of interference diagrams in the reduced-particle description used in our earlier work [6,9] diagram elements, shown outside the boxes in figs. 17 and 18, represent vertices for 2→3 and 2→ 4 particle transitions respectively. In the reduced-particle description of fig. 11b, these become effectively 0→1 and 0→2 transitions, with the states described by one transverse momentum or position variable (P or B) per effective particle.
The first expression for P |−iδH| in fig. 17 shows the general rule for converting ordinary matrix elements (the boxed formula) into the conventions of AI1 [6] for effectively 0→1 particle transitions. The only change here to the δH matrix element is that the normalization conventions of AI1 [6] include a factor of |2E n | −1/2 for each individual particle state. Our convention is that E n = x n E, where E is the energy of the original high-energy parent at the very start of the processes depicted in figs. 7-9, and where x n is the longitudinal momentum fraction of particle n with the convention that x n is negative for particles in the conjugate amplitude, so that n x n is always zero.
The final expression for P |−iδH| in fig. 17 shows what you get if you expressū 3 / ε * 2 u i in the same notation used in our previous work [6]. Here where the p's represent transverse momenta. For 3-particle (effectively 1-particle) states, the P ij are related by momentum conservation: 12 The P are related to the square roots of helicity-dependent DGLAP splitting functions. Here,

3)
12 For a discussion of (3.1), (3.2) and (3.4) in the context of our application and notation here, see specifically AI1 sections II.E and III [6].
FIG. 17: The box shows the Hamiltonian matrix element for e → γe, as in LCPT, normalized with relativistic normalization as in [13]. The next line shows the corresponding formula, in the normalization and conventions of AI1 [6], for an initial splitting e → γe in the amplitude (as opposed to conjugate amplitude), like the left-most (earliest time) vertex in fig. 11b or in figs. 7(a-c) and 8(h-k). As written, the latter formulas only apply to 2→3 particle transitions in the language of fig. 11a, which is equivalent to effectively 0→1 particle transitions in the language of fig. 11b. Note that these formulas show the matrix element of −i δH rather than δH, according to the convention of AI1 [6]. [The origin of the −i factor is the −i in the evolution operator e −iHt for amplitudes.] where z n ≡ x n /x i are the momentum fractions of the daughters relative to their immediate parent, P i→23 (1→z 2 , z 3 ) is the helicity-dependent DGLAP splitting function appropriate to the particle types and helicities of the initial and final particles in the splitting process, and the ± in the circular basis vector e x ±ie y is chosen accordingly based on those helicities. 13 The Fourier conjugate of P ij is, in our notation, The B|−i δH| formula in fig. 17 is simply the Fourier transform of the P |−i δH| formula. Now turn to fig. 18. The first expression for P 34 , P 12 |−i δH| there shows the general rule for converting ordinary matrix elements (the boxed formula) into the conventions of refs. [6,9] (AI1,ACI3) for instantaneous, effectively 0→2 particle transitions. In addition to the |2E n | −1/2 factors, there is an additional factor of |x 3 + x 4 | −1 that is associated with the normalization of the state |P 34 , P 12 that one obtains when reducing from a 4-particle description, of the state just after the interaction, to an effective 2-particle description. 14 There is nothing special here about the pairing of the indices 1234 in this normalization factor; one could just as well have written The only advantage to choosing P 34 , P 12 | in fig. 18 is that then the |x 3 +x 4 | −1 from the normalization factor there neatly combines with the (x 3 +x 4 ) −2 coming from the denominator fig. 18 is just the Fourier transform of the P 34 , P 12 |−i δH| formula.
We need not state whether our convention for light-cone components is The formulas for matrix elements, such as in fig. 18, give the same result either way.
From formulas such as figs. 17 and 18 for basic interactions, many variations follow from substitutions and complex conjugation. As an example, the interference-diagram rule in fig.  17 also gives the additional interference-diagram rules shown in fig. 19. This procedure has the advantage of definiteness, but one may also formulate a more generic rule that covers all the different variations-see appendix D.
To illustrate the possible relations, we've shown more variations in fig. 19 than we will actually need. The first line of vertex diagrams is not relevant for large N f , but at subleading order in 1/N f would appear in the virtual-correction diagram of fig. 20a. The second line of vertex diagrams would appear in the diagram for the complex conjugate of fig. 7a, shown in fig. 20b. But this diagram need not be evaluated separately since we will be taking 2 Re[· · · ] of the diagrams in fig. 7. The third and fourth lines of vertex diagrams in fig. 19 appear in figs. 8(i,k,m) and figs. 8(h,l), respectively. Note that these vertices correspond to (e → γe) * in the conjugate amplitude and not to eē → γ. They are drawn as they are just because the diagrams for the amplitude and conjugate amplitude have been sewn together. Depending on your taste, it might have been clearer to draw the third line of fig. 19 as fig. 20c, but we have chosen to keep with the same style of drawings as in ref. [6] (AI1).
Using these rules, and the additional vertex rules of appendix D, the calculation of the real double splitting diagrams of fig. 7 proceeds almost the same as our gluon-splitting calculations in our earlier work [6][7][8][9]. In particular, the calculation of the second line of diagrams in fig. 7, which involve the instantaneous 4-fermion vertex of LCPT, is very similar to the calculation of QCD diagrams involving the 4-gluon vertex [9]. Because these calculations are so close to previous work, we leave the details and analytic results to appendix E.
We'll mention just one qualitative detail here: During effectively 2-particle evolution (e.g. the middle shaded region of fig. 11b), the system evolves in medium like a coupled pair of 2-dimensional non-Hermitian harmonic oscillators, with two complex eigenfrequencies. 15 In (large-N c ) QCD double bremsstrahlung, both eigenfrequencies are non-zero. In QED, however, one vanishes. See appendix E for details, but this difference does not have any particular impact on the method of calculation.

B. An implicit approximation
There was an implicit approximation made in the above treatment of instantaneous LCPT interactions. Remember that the high-energy particles shown in our diagrams, such as fig (suppressing dependence on longitudinal momentum fractions, and working in the thickmedia limit being considered in this paper). In the limit of large E, this is parametrically small compared to any distance scale that characterizes the medium, such as the typical distance between medium particles. So, the chance that there is a medium interaction during the instant of the instantaneous interaction is parametrically small and can be ignored. 16 On a similar note, there are some other types of interactions in time-ordered LCPT calculations that we have also ignored. The LCPT rules that we have quoted assume that the LCPT Hamiltonian has been normal ordered. When normal ordering the instantaneous interactions, such as fig. 15, there are contractions that produce additional 2-point interactions in the normal-ordered Hamiltonian, which are often depicted diagrammatically [12] as in fig. 22. In our large-N f QED calculation, these can contribute to the virtual correction to e → γe as in fig. 23, which should be added to the diagrams we listed earlier in fig. 9. Fortunately, we can ignore these additional diagrams for reasons somewhat similar to those for the LCPT analysis of NLO deep inelastic scattering in refs. [19,21]. Because we are in the high-energy limit, we have ignored the masses of our high-energy particles. We will be using dimensional regularization, and normal-ordering contractions such as fig. 22 vanish in dimensional regularization for massless particles in vacuum, which is ultimately a consequence of dimensional analysis. Unlike refs. [19,21], however, our loops are not in vacuum, and the medium introduces a scale that in principle invalidates the argument that these loops will vanish. But because these particular loops involve an instantaneous interaction, and because transverse separations are suppressed by powers of the energy E, the effect of the medium on the loops of fig. 22 are negligible in the high-energy limit. So we may ignore the diagrams of fig. 23.

IV. VIRTUAL CORRECTIONS TO SINGLE SPLITTING
Consider the diagrams of fig. 8, which represent next-to-leading order corrections to single-splitting. Each of these diagrams contain a virtual eē loop. We shall now see that almost all of these diagrams are related in a relatively simple way to the double-splitting 16 In contrast, the effect of the medium on the evolution between vertex times in our diagrams (e.g. the shaded regions of fig. 11) cannot be ignored because the times between these vertices are of order formation times, which are parametrically large in the high energy limit.  fig. 7 (a,b,c,e,g respectively) by a simple diagrammatic procedure: Take the latest-time vertex in the diagram, and slide it around the right end of the diagram from being a vertex in the amplitude to being a vertex in the conjugate amplitude, or vice versa. For concreteness, a specific example is shown in fig. 24. We will call this a "back-end" transformation. Provided that the vertex being moved is the latest-time vertex in both diagrams, there is an extremely simple relationship between the values of the diagrams: Before integrating over any longitudinal momentum fractions, the two diagrams differ only by an overall minus sign. Heuristically, this minus sign can be roughly understood from the relation of virtual loops, through the optical theorem, to the probability of something not happening. So, for example, figs. 7(b,c) are, roughly speaking, related to the probability P 2 of one splitting e → γe later being followed by another, γ → eē. Figs. 8(i,j), on the other hand, are roughly related to the probability P 1 of one splitting e → γe not being later followed by another. 17 Any increase in P 2 should be accompanied by a decrease in P 1 by conservation of probability, and so one may expect these diagrams to be the same size but 17 This characterization makes sense for (1) a formal calculation of probabilities in the limit where the medium has some finite (but very large) size L, formally expanding to second-order in perturbation theory in hard splittings. That is, the limit α → 0 for fixed L and fixed particle energies (with L large compared to formation times and α here referring to the α at the scale characterizing the high-energy splitting). More precisely, this is the limit where the mean free time between splittings is large compared to L. But one may be interested in the opposite order of limits (2) L → ∞ with fixed α ≪ 1 and fixed particle energies-that is, the case where the mean free time between splittings is small compared to L. As discussed in detail in ACI2 [7], the formal calculation (1) can be used to figure out a correction to splitting probabilities that allows for the calculation (2), and so it is in the formal context of (1) that we always discuss our diagrams, even for those cases where our ultimate interest may be (2). with opposite signs. (In principle, the total probability that's actually conserved at this order is P 0 + P 1 + P 2 , where P 0 is the probability that neither emission occurs. But the cases of the diagrams that we are relating by back-end transformations are cases where changes to P 0 do not come into play.) By writing expressions for diagrams that are related by a back-end transformation, using the vertex rules of section III A and appendix D, and then tying them together with nparticle evolution in the medium as in fig. 11 and refs. [6,7] (AI1,ACI2), one may verify that back-end transformations really are that simple: there is just an overall sign difference. Note in particular that the longitudinal momentum fractions of the particles are the same in the two diagrams of fig. 24, for each of the three shaded regions of in-medium evolution. So those evolution factors match up identically between the two diagrams.
Virtual loops should be integrated over the longitudinal momenta of the particles in the loop. So the final relationship between the pairs of back-end related diagrams discussed above can be written as , where Γ is the rate, x e is the momentum of the daughter whose electron line is connected to the original electron, and y e is the other electron in the final state of real double splitting processes e → eēe, as in fig. 24. The integration limits come from the facts that (i) the longitudinal momentum fractions in the produced eē pair are y e and 1−x e −y e in these diagrams (relative to the initial electron), and (ii) longitudinal momentum fractions are always positive in LCPT.

B. Front-end transformations
There is a somewhat related relation between pairs of diagrams where we instead take the earliest-time vertex in a diagram and slide it around the left end of the diagram to move it from the amplitude to the conjugate amplitude or vice versa. We will refer to this as a front-end transformation, an example of which is shown in fig. 25. We'll discuss how to implement a front-end transformation in a moment, but first note its utility. Using a front-end transformation, along with complex conjugation as necessary, the γ → eē virtual correction diagrams of figs. 9(o,p,q,s,u) can be related to the real double bremsstrahlung diagrams of figs. 7(a,b,c,f,g), respectively. Using both a front-end and a back-end transformation on a single diagram, figs. 8(m) and 9(t) can be related to figs. 7(e,f), respectively, and fig. 9(r) can be related to fig. 8(k). In consequence, the only virtual diagram calculation that we will need to do from scratch is the boxed diagram of fig. 8(k): an (in-medium) photon self-energy loop in the "middle" of a single splitting e → γe.
Front-end transformations are more complicated to implement than back-end ones. If we were to simply slide the vertex around as in fig. 25 while keeping the labeling of all longitudinal momenta the same, we would get the transformation shown by the first two diagrams of fig. 26. There are two problems with the middle diagram in fig. 26, which is supposed to represent a virtual correction to γ → eē. The first problem is that our convention is always to let E refer to the energy of the original high-energy particle in any process, which in this case would be the photon. The second is that x e is positive in the diagram of fig. 26, which is then inconsistent with the fact that the longitudinal momentum x e E y e E vs. x e E shown in the second diagram should be negative, given our convention that longitudinal momenta of conjugate-amplitude particles are negative in our interference diagrams. This second problem can also be visualized by comparing the earlier figs. 24 and 25 for back-end and front-end transformations. In the back-end case of fig. 24, both the original diagram and the transformed diagram had exactly the same particles evolving in each of the three shaded regions. In the front-end case given by the first two diagrams of fig. 25, the last shaded region is the same, but an "amplitude" particle (blue line) has been switched to a "conjugate amplitude" particle (red line) in each of the first two shaded areas. That means that the evolution in those shaded areas is not exactly the same for those two diagrams, and so their relation is not as simple as it was for a back-end transformation.
There is a simple way to overcome both obstacles, which is to make the change of variables when making a front-end transformation like fig. 25. This converts the momenta of the middle diagram of fig. 26 into those of the last diagram. It also negates the value of x e . As long as one writes formulas for evolution in diagrams in a way that is general enough to correctly handle the cases of both positive or negative longitudinal momentum fractions, similar to ref. [6] (AI1), then all will be well. The transformation (4.2) is its own inverse, and so works just as well transforming γ → eē back again to e → eēe. In terms of the diagram labels in figs. 7 and 9, is the front-end analog of (4.1). 18 Here we have taken 2 Re[· · · ] of all the diagrams, which we must do anyway at the end of the calculation, and which here obviates the need to specify exactly which diagrams need to be complex conjugated after the front-end transformation. note that we have written the γ→eē rate as dΓ/dy e to conform with our convention that the electron produced by pair production is labeled y e , as in the last diagram of fig. 26. Performing both a front-end and back-end transformation relates 2 Re dΓ dy e (t) = +2 Re As shown in fig. 27, the additional step of (x e , y e ) → (1−y e , x e ) is needed to make the labeling of the momentum fractions in diagram (t) match up with our conventions. A different variant of the front-end transformation is needed for interference diagrams whose earliest vertex is an instantaneous 4-fermion interaction, such as fig. 7(e). Fig. 28 shows diagrammatically the difference between (i) "simply sliding the vertex around" while keeping the labeling the same vs. (ii) the actual result we want to achieve with a front-end transformation. We can implement the necessary relabeling of momenta by making the change of variables 5) vs. C. The fundamental virtual diagram

Overview
As previously discussed, the one virtual correction diagram which cannot be related to the e→eēe diagrams is fig. 29, which is a more detailed version of fig. 8k. Following notation similar to that of ref. [6] (AI1), we will refer to this as the xyyx diagram since it involves, in order, (i) emission of an x e electron in the amplitude, (ii) emission of a y e electron in the amplitude, (iii) re-absorption of the y e electron in the amplitude, followed finally by (iv) emission of the x e electron in the conjugate amplitude. This virtual correction diagram contributes dI dx e xyyx = 1−xe 0 dy e dI dx e dy e xyyx (4.7) to the differential emission probability for e → γe, where, adapting the notation of refs. [6,7] (AI1,ACI2), In broad outline, this formula simply convolves the relevant vertex matrix elements at different vertex times with medium-averaged propagators between those times. 19 The notation [right]. Our notation for vertex times is also shown. (3.4); we use the the letter C just to help distinguish B ij variables used in the context of 4-particle (effectively 2-particle) evolution. The indices on C ij refer to the particles involved in the 4-particle evolution according to the convention of the right-hand drawing in fig. 29, for which we introduce the labeling [6] ( of longitudinal momentum fractions. Initially, the evaluation of (4.8) proceeds very similarly to that of other diagrams, following refs. [6,7] (AI1,ACI2). One puts in explicit formulas for the δH matrix elements and for the time-evolution factors. It is then possible to analytically perform all the integrations implicit in (4.8) over the intermediate transverse positions specified by the B's and C's, as well as all but one of the time integrations. Due to time translation invariance, one combination of those time integrations just gives a factor of total time T. After all these integrations, the differential rate dΓ ≡ dI/T corresponding to (4.8) may be reduced to a single remaining integral over the time difference ∆t between the two intermediate vertex times in the process under consideration. Schematically, For the xyyx process of fig. 29, ∆t ≡ t y ′ −t y is the duration of the photon self-energy loop. The integrand F (∆t) has a somewhat complicated formula (schematically similar to other diagrams), which we leave to appendix F. Here, in the main text, our goal is to give a broad overview of what is involved in dealing with the UV divergence of the ∆t integral in (4.10). To be concrete: For small ∆t, the integrand turns out (see appendix F 2) to have non-integrable divergences where the magnitude of can be interpreted physically as the scale of the QED inverse formation time (2.4) for leadingorder e → γe. UV divergences from ∆t→0 in the calculation of time-ordered interference diagrams are not restricted just to loop diagrams like fig. 29. Individually, the real double splitting diagrams of fig. 7(a-c) for e → eēe also have ∆t→0 divergences when written in the form d(∆t) F (∆t), analogous to the detailed discussion of the QCD case in refs. [6][7][8] (AI1-ACI3). Those divergences cancel between the three diagrams of fig. 7(a-c), also similar to refs. [6][7][8]. Consistently evaluating the remaining finite part is a little tricky and was accomplished [8] by computing the individual diagrams using dimensional regularization.
Because the virtual diagrams of fig. 8(h-j) are related to fig. 7(a-c) by a simple back-end transformation (4.1), the divergences of (h-j) will also cancel. Additionally, we find that diagrams involving instantaneous longitudinal photon interactions are all individually finite. That means that the only uncanceled UV divergence in our calculations of e → eēe and e → γe is the one for the xyyx diagram shown in both figs. 7(k) and 29. [The divergence of the γ → eē diagram of fig. 9(r) is similarly uncanceled, but, as already noted, that diagram is related to the other by front-and back-end transformations.] Unfortunately, our previous method [8] for calculating diagrams in dimensional regularization breaks down for the xyyx diagram, for reasons that will be explained later. We will leave the (we think interesting) details to appendices. Here we convey the gist of the method.

Dimensional Regularization: Strategy
In this paper, the symbol d ≡ d ⊥ will represent the number of transverse dimensions, not the total number 4−ǫ of space-time dimensions. It might reduce the opportunity for confusion if we always wrote our d as d ⊥ , but we drop the subscript to save space in what will be some complicated equations. We do not know how to do the entire integral (4.10) for arbitrary d. Fortunately, we do not have to. The only divergence comes from ∆t → 0, and so we only need to regulate the small-∆t portion of the integral. Let a be an arbitrarily small but finite time scale, and split (4.10) up into (4.14) We only need dimensional regularization for the first term on the right-hand side. Specifically, where F d (∆t) is the integrand for d transverse dimensions. Roughly speaking, our goal is to use small-∆t approximations to do the ∆t < a integral analytically, and then use numerical integration to do the ∆t > a integral. But there is a problem: As we take a → 0 (in order to make our calculation of the first integral arbitrarily accurate), the second integral ∞ a F 2 (∆t) will blow up, which is undesirable for numerical integration. We can isolate this problem by rewriting the right-hand side of (4.15) as where D 2 (∆t) is any convenient function that • matches the divergence of F 2 (∆t) as ∆t → 0 [which is proportional to (4.11)]; • 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 (4.16) is then convergent by itself in the a → 0 limit, and so we may replace (4.16) by Our strategy will be to evaluate the first two integrals in (4.17) analytically in the limit of small a (and verify cancellation of the a dependence) and to evaluate the (convergent) third integral numerically. The choice of D 2 (∆t) is not unique, but the final answer for (4.17) will be independent of the details of that choice. The particular choice we will find convenient is which reduces to (4.11) for small ∆t but also falls off quickly enough as ∆t → ∞ [because csc(Ω ∆t) falls exponentially for the Ω ∝ √ −i of (4.12)]. But there is no difficulty if the reader thinks some alternative choice would have been simpler or more natural.
We placed scare quotes around the "a→0" limit in (4.17) because we should clarify how it coordinates with the ǫ → 0 limit of dimensional regularization. The earlier split (4.15) of the ∆t integration is invalid unless we can set d=2 for ∆t ≥ a, up to corrections that disappear as ǫ → 0. One feature of dimensional regularization is that it converts integer power-law divergences into non-integer power-law divergences. For example, d(∆t)/∆t might become proportional to d(∆t)/(∆t) 1−ǫ or d(∆t)/(∆t) 1−2ǫ . (We will see details later.) When integrated over ∞ a d(∆t), the results will correspondingly have non-integer power-law dependence on the lower cut-off a, and we will encounter factors like (aΩ) ǫ , where Ω is the relevant frequency scale of the problem. We need to be able to expand such a factor in small ǫ as 1 + ǫ ln(aΩ) + · · · in order to recover the d=2 result for the ∆t ≥ a integral, which means that ǫ ln(aΩ) will need to be small. The "a→0" limit in (4.17) is therefore notational short-hand for expanding expressions assuming that the small a has been chosen relative to the (eventually) arbitrarily small ǫ such that Before moving on to discuss the application of these techniques to the xyyx diagram of fig. 29, we should clarify how we plan to handle photon and fermion helicities in dimensional regularization. We will implement the Conventional Dimensional Regularization (CDR) scheme with MS renormalization. These are, for instance, the conventions used by the Particle Data Group [23] and the community at large to quote, compare, and take worldaverages of determinations of α s (M Z ). For a nicely brief summary of different flavors of "dimensional regularization" used by some authors, see, for example, section II.A of ref. [22]. 20

Dimensional Regularization: Application to xyyx
The δH matrix element in (4.8), which are the vertex factors of figs. 17 and 19 and appendix D, are all proportional to the transverse momentum P characterizing the splitting, which is equal to the P ij of any pair of the particles directly involved in the splitting (or merging). Below, take this to be the P ij of the two daughters. The proportionality of δH matrix elements to P holds in any dimension d. The proportionality constant, which involves DGLAP splitting amplitudes and consideration of particle helicities, depends on d, but for the moment we won't keep track of those details. In transverse-position space, the factors of P become factors of ±i∇δ (d) (B). Plugging this into the starting formula (4.8) for the xyyx diagram, and using the δ functions to perform the implicit integrations over the B variables, gives where the various gradients ∇ are contracted in ways that depend on the helicity-dependent DGLAP splitting functions not explicitly shown here. The integral over the relative vertex times in the xyyx time-ordered diagram of fig. 29 can be written as In the multiple scattering (q) approximation, the 3-particle (effectively 1-particle) propagators B ′ , t ′ |B, t are equivalent to those of a single d-dimensional harmonic oscillator in non-relativistic quantum mechanics, whose mass M and (complex) frequency Ω are related 20 For the quantities computed in this paper, there will be no difference between using what ref. [22] identifies as CDR and HV dimensional regularization, once we write final formulas for differential rates in terms of MS-bar renormalized α and set ǫ = 0. In particular, for the quantities we calculate, we will not have any IR/collinear divergences that need to be canceled between real and virtual emission diagrams, and so it does not matter if we treat real and virtual particles on exactly the same footing in d=2 − ǫ dimensions. In part this is thanks to the fact that medium effects cut off infrared and collinear singularities for splittings contained within the medium.
to the longitudinal momentum fractions of the three particles. In the particular case of the xyyx diagram, the momentum fractions are the same for both regions of 3-particle evolution and give the Ω i of (4.12) and (4.22) The formula for a harmonic oscillator propagator, is simple enough that the integrals over the initial and final times in (4.20) can be carried out explicitly [8], giving 21 where now it is the two Bs and two ∇s that are contracted in ways that depend on the helicity-dependent DGLAP splitting functions. The 4-particle (effectively 2-particle) propagator C ′ 41 , C ′ 23 , t ′ |C 41 , C 23 , t is also a harmonic oscillator propagator, but this time for a coupled set of two harmonic oscillators. A number of terms are generated when one takes derivatives ∇ C ′ 23 ∇ C 23 of this propagator as in (4.24), but only one of them leads to divergences of the ∆t integral when ǫ = 0, and that is then the only term we need to treat with dimensional regularization. Details are given in appendix F. The small-∆t behavior of that term, which is all we need to regulate to implement (4.17), is specifically where the exponential factor is analogous to the exponential factor in the single harmonic oscillator propagator (4.23) and corresponds to the exponential of the double harmonic oscillator problem after setting C ′ 23 = 0 = C 23 . The small-∆t expansions of the coefficients X in the exponent turn out to be [The fact that the expansions to this order can be written in terms of Ω i (4.12), without reference to the actual eigenfrequencies Ω ± of the double harmonic oscillator problem, is non-trivial. See appendix F 2 b.] In (4.25), we also now show explicitly the particular combination of DGLAP splitting functions that appear, which for this divergent term happen to combine into a simple product of the (d=2−ǫ dimensional) spin-averaged splitting functions P γ→e for e → γe and γ → eē.
The technical problem we face to carry out dimensional regularization is how to do the integrals in (4.25) with (4.26) in the small-a limit (4.19). For similar integrals in ref. [8] (ACI3), we were able to argue that the exponential factor in the integrand limited the range of the Bs enough in the small-∆t limit that one could make small-argument expansions of the Bessel functions K d/4 ( 1 2 |M i |Ω i B 2 ). Unfortunately, that is not the case here. If the exponent in (4.25) were evaluated with (4.26) at leading order in small ∆t, the exponential would be exp  The integrand remains unsuppressed when the arguments of the Bessel function are of order 1, and so small-argument expansions of the Bessel functions are not applicable. Having to deal with the Bessel functions as they are makes the integral harder to do.
In appendix G, we show how to evaluate the integral which is proportional to the one shown in (4.25), to obtain 31) where γ E is Euler's constant. The a dependence will be canceled, as it must, when one adds together all the terms of (4.17) for the calculation of xyyx. When we later combine the xyyx diagram with leading-order e → γe, the 1/ǫ pole will be absorbed, as it must, by the known renormalization of α EM in QED.
The rest of the calculation of xyyx is mostly just a matter of combining the above with calculations and formulas that are closely analogous to those for diagrams analyzed in our previous work [6][7][8]. We will write our decomposition (4.17) of the calculation as

Renormalization
There are two equivalent ways to think about renormalization and the ultimate cancellation of the xyyx UV divergence appearing in (4.31). One is renormalized perturbation theory: Add counter-term diagrams such as fig. 30, so that the combination dΓ dx e xyyx + dΓ dx e counterterm diagram (4.33) is finite. Note that in the large-N f limit of QED, there is no electron wave function renormalization and no vertex renormalization, and in QED the photon wave function renormalization is equivalent to charge renormalization. The other viewpoint, which we find more straightforward to implement, 22 is to imagine that the contribution to dΓ/dx e from all e→γe diagrams, including the leading-order process represented by the xx diagram of fig. 31 (plus its complex conjugate), are initially computed 22 The second viewpoint avoids having to sort out a few possible sources of confusion. For example, when computing amplitudes in quantum field theory, one multiplies external legs by Z 1/2 field instead of Z field , but it is the latter that comes from the divergences of self-energy loops such as the photon self-energy loop in the xyyx diagram. The seeming difference can be resolved by realizing that the leading-order xx diagram represents a rate rather than an amplitude, and so should be multiplied by |Z 1/2 field | 2 = Z field , which generates the counter-term diagram in fig. 30. But one could also worry about why we have not included other time-orderings of that counter-term diagram, such as ones where the counter-term occurs later than both emission vertices. Those contributions turn out to cancel between diagrams with the counter-term in the amplitude and diagrams with the counter-term in the conjugate amplitude, because of the different signs in the corresponding evolution operators exp(∓i δH t). The second viewpoint on renormalization allows us to bypass all of these considerations. using the bare coupling instead of the renormalized coupling. Afterwards, we convert to the MS-renormalized coupling using the known relation which is equivalent to fig. 30. Note that, because it is multiplied by 2/ǫ, we will need to use a d=2−ǫ formula for the leading-order [dΓ/dx e ] xx above, expanded through O(ǫ). The explicit formula is given in (F45) of appendix F. Now combine the decomposition (4.32) of xyyx into calculable divergent and finite parts with the renormalization (4.37), and choose D 2 (∆t) according to (4.18). We find (see ap- where µ is the renormalization scale. [dΓ/dx e ] (subtracted) is as defined earlier in (4.32), with the result given by (F41) of appendix F. An aspect of handling xyyx divergences worth mentioning here, that was unnecessary in our previous application of dimensional regularization to real double splitting g → ggg in ref. [8] (ACI3), is that we find we need the d = 2−ǫ result for a DGLAP splitting functionspecifically P (d) γ→eē , which is the additional splitting function factor in the xyyx divergence (4.25) compared to the leading-order xx diagram of fig. 31. where [dΓ/dx e dy e ] (k) represents the y e -integrand that integrates to (4.38), which may be interpreted as dΓ dx e dy e (ren) with the understanding that 1 0 dy e δ(y e ) = 1. Above,Ω i ≡ Ω i sgn(M i ), which is needed to generalize the xyyx result to work with front-end transformations, as explained in appendix F 2 c.
Since we combine leading and next-to-leading order results when we renormalize as in (4.36), we note here that the relation between leading-order results for γ → eē and e → γe

V. CONCLUSION
The complete list of final analytic results for diagrams is sprinkled throughout this paper but consists of everything in appendices E 1 and E 2, plus eqs. The work in this paper completes the analytic calculation of the effect of overlapping formation times in QED in the simplifying parametric limit α ≪ N f α ≪ 1. Even though that is a QED limit only a theorist could love, we think it will be interesting to look at the size of those effects on shower development as one changes the value of N f α to be closer and closer to O(1), as a first exploration of how to carry out similar calculations for QCD if one pushes N c α s (Q ⊥ ) to be larger and larger to answer the question of how small or overwhelming overlap effects are in practice. We leave that QED analysis, as well as discussion of numerical results for the formulas presented here, to forthcoming work [15]. With these QED results in hand, it should now be possible to complete similar calculations for large-N c QCD.
Refs. [20,22] study processes somewhat similar to ours but in the context of next-toleading-order DIS in the dipole approximation appropriate to studying small-x physics. As is standard, one can use the optical theorem to relate DIS cross-sections to the self-energy of the virtual photon as depicted at leading order, for example, by the time-ordered diagram in fig. 32a. In their application to small-x physics, the medium is very thin compared to the formation time, and we've depicted its extent by the very thin gray region in fig. 32a. Fig. 32b shows the same process but now drawn using the conventions that we have used in fig. 7. 24 In this paper, however, our explicit calculations are for the case where the medium is thick compared to the relevant formation times, and so would be analogous to fig. 32c  rather than fig. 32b.
The photon in DIS is virtual, with a virtuality Q 2 that should not be ignored. In contrast, in our calculations of in-medium showering, we approximate the initial high-energy particle 24 In case the reader is wondering how the overall sign matches up between fig. 32a and 32b: In our formalism, the red portion of the diagram represent the conjugate amplitude and so evolves with e +iHt instead of e −iHt . When doing time-ordered perturbation theory, the sign difference in e ∓iHt manifests as a sign difference between red vertices and blue vertices in our diagram. in figs. 7-9 as on-shell, with negligible mass or virtuality.
Yet another difference is that refs. [20,22] write their final answers explicitly in terms of Wilson lines running through the medium, which should be averaged over medium fluctuations. Because the medium is thin, these Wilson lines do not have time to move transversely as they cross the medium. For thick media, the paths inside the medium do move transversely, and this dynamics is incorporated in our treatment of the LPM effect [6] (based on Zakharov's picture [18]) in the language of two-dimensional quantum mechanics with a non-Hermitian potential energy. This is an approximation that relies on correlation lengths in the medium being small compared to both formation times and the thickness of the medium. Finally, our explicit calculation in this paper further takes this potential to be given by thê q approximation, appropriate for the LPM effect at high energy (but see appendix C).

Appendix B: Leading-order splitting formulas
In this appendix we quote the leading-order LPM results for e → γe and γ → eē for reference. As in the rest of this paper, we take the high energy limit, assume a thick medium, and make the multiple scattering (q) approximation. As discussed in appendix C, these results, in different notation, go all the way back to Migdal [2].
1. e → γe In our notation, the leading-order e→γe rate is 25 and In terms of the xx diagram of fig. 31, the above result is 25 QCD versions go back to BDMPS [24,25] and Zakharov [18]. To relate the QED and QCD results in our own notation: (B1) is AI (1.5a) with P g→g (x) replaced by P e→e (x e ). The formulas for Ω 0 in QED and QCD are both special cases of the general formula Ω 2 = −i[(q 1 /x 1 ) + (q 2 /x 2 ) + (q 3 /x 3 )]/(2E) of AI (2.33b), whereq i is zero for the photon.

γ → eē
Leading-order pair production γ → eē is essentially the same except that and Appendix C:q in QCD and QED In this appendix, we mention some qualitative differences between QCD and QED concerningq and its logarithmic dependence on energy at high energy. Though the development of this paper does not depend on the medium itself being weakly-coupled, we start by discussing that case.

1.q for weakly-coupled systems
For weak coupling,q is given byq where dΓ el /d 2 q ⊥ is the rate of elastic scattering from the medium for momentum transfer q ⊥ perpendicular to the direction of motion of the high-energy particle. Coulomb interactions give for q ⊥ large compared to the inverse electric screening length ξ −1 of the medium. 27 For fixed small coupling α, this behavior leads to a logarithmic UV divergence in the integral (C1) forq in both QED and QCD:q One may then defineq(Λ) as a function of some UV cut-off Λ on q ⊥ . It turns out that the LPM calculation of the rate for splitting processes such as hard bremsstrahlung depends onq(Q ⊥ ) where Q ⊥ is the order of magnitude of the total q ⊥ transferred to the high-energy particle during a formation time. As we'll review below, this effect (though not in this language) goes all the way back to the results of Migdal [2]. 28 That fact notwithstanding, it is interesting that QCDq(∞) as calculated by (C3) is actually convergent if one accounts for the running of the coupling as α = α(q ⊥ ) in this particular calculation. The slow decrease of α with increasing momentum due to asymptotic freedom is just enough to then make the integral (C3) convergent. 29 In contrast, the running of the coupling in QED makes the divergence slightly worse since α EM (q ⊥ ) grows with momentum.
But this apparent conclusion thatq(∞) is finite for QCD, however, is an artifact of ignoring other higher-order corrections (besides running of the coupling) to the leadingorder analysis (C3), as we now review.

Small-x logs in QCD
In this paper, we have looked at splitting processes such as high-energy bremsstrahlung, which are suppressed by a factor of α s (Q ⊥ ), which is small at high enough energy even if the medium itself is strongly coupled, i.e. even if α s (ξ −1 ) is large. Liou, Mueller, and Wu [30] have shown that the contribution of high-energy gluon bremsstrahlung toq in QCD is enhanced by a double logarithm, arising from emission of nearly-collinear gluons with longitudinal momenta between the medium scale and the high-energy scale E. This double log is related to double logs that occur in small-x physics. At large energy, the double logarithm compensates for the small α s (Q ⊥ ) and so QCDq(∞) is divergent after all, even after resummation of leading logarithms at all orders. Refs. [3][4][5] extended this analysis to the use ofq in the context of calculations of high-energy splitting rates. Ref. [3] argues that the effective value ofq in splitting calculations scales like L 2 √ᾱ s for large energy E, with α s ≡ N c α s /π, and where L in the present context (infinite medium) means the formation length. Since L scales like L ∼ E 1/2 (for fixed x and smallᾱ s ), that givesq eff ∼ E √ᾱ s , which diverges as E → ∞.

The upshot for this paper
Whether we are talking about QED or QCD, there is some sort of logarithmic dependence ofq(Q ⊥ ) on energy. Theq in this paper should be fixed to the one appropriate for the energy E. By making the approximation that we can describe scattering from the medium in terms of a fixed effective value ofq, we are ignoring one class of sub-leading corrections. In the long term, one should be able to do a more complete calculation. In the medium term, we will sidestep this particular issue in forthcoming work [15] by looking at certain characteristics of high-energy showers that are not sensitive to the precise value ofq.
Migdal [2] for high-energy QED showering off of a medium made up of atoms, a useful summary of which can be found in the review by Klein [31]. Migdal's explicit solutions to his equations rely on assuming that Coulomb logarithms are large, and he does not try to precisely compute the constants under the logarithms. In the case of significant LPM suppression (what Migdal would call s ≪ 1), Migdal's results for e→γe can be rewritten in the following form at the same leading-log order: 30 with P e→γ (x γ ) the DGLAP splitting function Physically, the argument ofq eff in (C5) represents Q ⊥ ∼ √q t form (from the definition ofq), remembering that the frequency Ω is of order 1/t form and that rough approximations are all that are needed for the argument of the logarithm in (C6) for a leading-log analysis. Theq in the argument ofq eff can be interpreted as self-consistentlyq eff itself, but it does not matter at leading-log order. In (C6), our notationq eff (b −1 ) is motivated by b ∼ 1/Q ⊥ because we find it more convenient to express the right-hand side of (C6) in terms of transverse distance scales. Above, m e is the electron mass, R A is the nuclear radius (which Migdal somewhat obscures by approximating R A ≃ 0.5Z 1/3 α/m e ), a Z is the length scale for screening of the nucleus's Coulomb field by atomic electrons (for which Migdal uses the Thomas-Fermi approximation a Z ≃ Z −1/3 a 0 ), and Z is atomic number. The case b ≫ m −1 e not shown above would correspond to no LPM suppression (what Migdal would call s ≫ 1), and the case b ≃ m −1 e would correspond to the transition where the LPM effect is first turning on (which Migdal would call s ∼ 1).
The point of writing these formulas in the above form is that (C4) is precisely the leadingorder result (B1) in our notation. And the leading-log calculation (C1) ofq in Migdal's application iŝ where n is the density of atoms. There are now two cases to consider. (i) We mentioned previously that the scale Q ⊥ acts as a UV cut-off on the relevant value ofq eff for splitting 30 It is easiest to take Migdal's results from Klein (72-77) [31]  calculations. Taking Λ ∼ Q ⊥ , and then rewriting Q ⊥ as 1/b for the sake of expressing scales in terms of transverse distance scales instead of transverse momentum, gives the logarithm shown for the R A b ≪ m −1 e case of (C6). (ii) Ignoring the sub-structure of the nucleus (whose effects are suppressed by powers of Z −1 except at very much higher energies 31 ), there is a UV-cutoff on the nucleus's Coulomb field at b ∼ R A . The effective cut-off for splitting calculations is therefore Λ −1 ∼ min(Q ⊥ , R −1 A ) instead of just Λ ∼ Q ⊥ , which accounts for the other case of (C6).

Appendix D: Diagrammatic vertex rules
In the main text, figs. 17-19 gave examples of LCPT vertex rules and the corresponding rules in our formalism. The rest of the rules that we need for the large-N f QED diagrams of figs. 7-9 are shown explicitly in figs. 33-36. In some transformations, we must negate momentum variables (x n , p n ) → (−x n , −p n ), as indicated, to account for our convention that momentum variables are negated for particles in conjugate-amplitudes (red lines in the diagrams). Transverse position variables b n are unaffected. When applied to the definitions (3.1) and (3.4) of P ij and B ij , this transformation takes (P ij , B ij ) → (P ij , −B ij ).
The annoying part of working with the rules laid out for matrix elements is that one must be careful about Fermi statistics in the representation of states, e.g. |eē = −|ēe . The boxed LCPT rule forēe → γ in fig. 34, for example, could just as well have been written as f|δH|3, 2 rel = +gv 2 / ε f u 3 . The sign would be compensated, in a photon self-energy diagram for example, by whether the propagator between the two vertices was 2, 3; t ′ |2, 3; t or its negative 2, 3; t ′ |3, 2; t . (The minus sign in the last case would be equivalent to the usual accounting where one says that fermion loops come with minus signs.) We should also clarify a point about our formulas for e → eēe matrix elements, such as 2, 3, 4|δH|i rel and related formulas in figs. 18 and 35. In these cases, we are only giving the specific contribution associated with the accompanying vertex diagram. For example, the contribution from swapping the two final-state electrons (which in any case is sub-leading in 1/N f ) is not included, and would correspond to drawing a different vertex diagram where the "2" and "4" lines were switched.
Keeping track of all the signs and different cases for matrix elements is painstaking. However, there is an equivalent way to formulate our rules for diagrams that makes it easier and a little bit more like the conventions for Feynman rules. All the cases of 1↔2 splittings can be subsumed by the rule shown in fig. 37, supplemented by a minus sign for each fermion loop in the interference diagrams of figs. 7-9 and 31. The rule of fig. 37 sometimes differs by a sign from the particular conventions of the −i δH and +i δH matrix elements we have written down for individual vertices, but the supplemental minus signs for each fermion where R p is the proton radius, the charges inside the nucleus do not contribute coherently to scattering. In this case, the contribution from R −1 A q ⊥ b −1 to the q ⊥ integral (C7) forq eff is at most of order 8πZαn ln(R A /R p ) (where Zn is the density of protons in the system), which is suppressed by a power of Z and so small (for moderate to large Z) compared to the contribution from a −1 Z q ⊥ R −1 A shown in the first case of (C6). If one goes to energies high enough to probe the substructure of the nucleons (which Migdal did not know about), then, as a matter of principle, eventually the contribution of scattering from individual quarks would become important at sufficiently small b ≪ R p .    FIG. 34: Like fig. 33 but for inverse pair production.    loop brings the two different procedures into agreement and also would produce the correct Fermi-statistics sign for exchange diagrams such as fig. 38 (which are sub-leading in 1/N f ). There are different ways one could assign signs to the vertices that would give the same net overall sign for interference diagrams such as figs. 8-9; we've picked one of them. The rule of fig. 37 has the same form as the similar rule given in AI1 [6] for the 3-gluon vertex except for details about signs. In the rule, the factor Pē eγ (x i , x j , x k ) depends implicitly on the helicities h of the lines as measured in the directions of the small arrows. The definition of Pē eγ is as in AI1 section IV.E, adapted here for the QED case as with spin-dependent DGLAP splitting functions and e ± ≡ e x ± e y . The zeros above are a consequence of chirality conservation. As in AI1 [6], the P (x i → x j x k ) are defined in terms of the usual DGLAP splitting functions by where z ≡ x/x i are the momentum fractions of the daughters relative to their immediate parent. The advantage of the P (x i , x j , x k ) is that they are normalized so that they are k i j symmetric with respect to permuting the parent with the daughters and so are the same for e → γe and γ → eē, as indicated in (D2). Eqs. (D2) only show half of the helicity cases; the other half are given by P h 1 ,h 2 ,h 3 = P −h 1 ,−h 2 ,−h 3 and so Similar rules for the instantaneous photon interactions are given by fig. 39, which is similar in implementation to the 4-gluon vertex rule of ACI4 fig. 10 [9]. The rule here is normalized to cover all the cases in this paper for 0 ↔ 4 or 2 ↔ 4 particle transitions, which involve instantaneous e → eēe in either the amplitude or conjugate amplitude.
[We have not bothered presenting a generalized rule for 3 → 3 transitions such as in figs. 7(g) and 8(n) for e → eēe and 9(u) for eē → eē, since these diagrams give zero by symmetry arguments similar to those of ACI4 section III.B [9].] will be ∆ dΓ dx e dy e = ∆ dΓ dx e dy e seq + dΓ dx e dy e (I) + dΓ dx e dy e (II) , where the pieces are given below and "I" indicates a contribution involving an instantaneous 4-fermion interaction.
In large N f , we distinguish the final-state electron that carries the flavor of the initial-state electron and refer to it as x e . The electron in the pair produced from the photon (which has a different flavor in large N f ) is referred to as y e . Because we distinguish these particles in our large-N f formulas, ∆dΓ/dx e dy e is normalized so that, in applications, final-state integrations should be performed as without any factor of 1 2 for identical final-state particles.

Sequential diagrams a. Generic Formulas
The calculation of the diagrams of fig. 7(a-c) mimics that of the corresponding diagrams in ACI2 [7]. Using the same notation, the results can be directly taken over, except for some normalization factors and low-level details that we will get to in a moment. Adapting the results summarized in ACI4 appendix D.3 [9] 32 , the sum of figs. 7(a-c) is ∆ dΓ dx e dy e seq e→eēe = 2N f A seq (x e , y e ) (E3) 32 The summary in ACI4 appendix D [9] incorporates some corrections [8] to the analysis of "pole" terms in the earlier paper ACI2 [7]. So it is better to take formulas from ACI4 appendix D than directly from the original ACI2 analysis.
with A seq (x e , y e ) ≡ A pole seq (x e , y e ) + +∞ 0 d(∆t) 2 Re B seq (x e , y e , ∆t) + F seq (x e , y e , ∆t) , (E4) A pole seq (x, y) = and There are a few small differences so far with the g → ggg formulas.
• The naive translation of C 2 A α s to QED is , which is 1 in QED, and (ii) there are N f possible flavors of the eē pair produced by the virtual photon in e → γ * e → eēe.
• In the QED analysis, we do not have to worry about large-N c "color routings," which affected the choice of how to write the g → ggg version of (E3), as described in ACI2 section 2.B.1 [7]. Each QCD diagram was the sum of two color routings, with each routing represented by A. If we adopt the translation (E11), that means that A will represent half of each QED diagram, which is the reason for the overall factor of 2 in (E3). We could have instead absorbed this factor of 2 into the normalizations of D seq , F seq and A pole seq , but it seemed more convenient to normalize the lower-level formulas in exactly the same way as the QCD case.
• For g → ggg, there are three identical particles in the final state, and refs. [6,7,9] correspondingly add together all permutations of A(x, y) corresponding to permuting x, y, and z ≡ 1−x−y. For e → eēe, there are two identical particles in the final state in ordinary N f =1 QED, and so one might reasonably think that, analogously, the A seq (x e , y e ) on the right-hand side of (E3) above should be A seq (x e , y e ) + A seq (y e , x e ). If so, one would integrate over final-state particles as 1 2 1 0 dx e dy e θ(1−x e −y e ) in applications of dΓ/dx e dy e , where the factor of 1 2 would avoid double counting of identical final states. In large N f , however, we distinguish the two electrons as explained previously. To allow applications the option of tracking the fate of the initial-flavored electron, we have chosen not to include the x e ↔y e permutation in (E3). When using our large-N f dΓ/dx e dy e in applications, one should correspondingly integrate as in (E2), without any final-state factor of 1 2 . The quantities I seq n in (E7) are defined the same way as in ACI4 appendix D.2 [9], which we repeat here for reference.
The formulas for (X, Y, Z) seq in (E7) in terms of the eigenfrequencies and eigenmodes of the problem are also given in ACI4 appendix D.2 [9], but we will see below that one of those eigenfrequencies, Ω − , vanishes in our QED application. The formulas then specialize to X seq where is a matrix of appropriately normalized modes of the double harmonic oscillator problem in the basis (C 41 , C 23 ) used by ACI2 [7] at the y vertex of the xyxȳ diagram, and is a permutation appropriate to thex vertex.
b. QED formulas for (ᾱ,β,γ) As in ACI2 [7], the (ᾱ,β,γ) in (E5) are functions of x e and y e that represent various combinations of helicity-dependent DGLAP splitting functions from the vertices of fig. 7a. The relevant QED splitting functions are different from those of the g → ggg process in QCD. Performing the same calculation as in ACI2 appendix E [7], but using the QED splitting functions for e → γe and γ → eē appropriate to fig. 7a, we find 33 where Note thatγ = −ᾱ.
One may check thatᾱ where which is the straight-forward translation to QED of a similar relation for g → ggg. 34

Frequencies Ω and eigenmodes
For all of the e → eēe diagrams of fig. 7, we will need the relevant frequencies Ω for 3-particle and 4-particle evolution, and the eigenmodes for 4-particle evolution. 33 Specifically, this is the e → eēe analog of ACI2 eq. (E4) [7]. The use of the letter z here is unrelated to the use in (D3). The absolute value signs on |x e y e z e | may seem redundant here, but they are included for a reason similar to footnote 21: to make sure that (ᾱ,β,γ) behave appropriately under front-end transformations (4.2). With the absolute value signs, a front-end transformation maps (ᾱ,β,γ) into (1 − x e ) 10 times the analogous helicity-averaged product of P's that one would have constructed for the last diagram of fig. 26. 34 See ACI2 eq. (E5) [7]. This is a re-assuring check because it has to hold in order for sequential diagrams to match up with sequential "Monte Carlo" when the two splittings are far separated in time. That requirement is implied by ACI2 footnote 28 [7] and the need for ACI2 eqs. (C7) and (C13) to match up accordingly. Analogous statements must hold for QED.

a. 3-particle evolution frequency
Quite generally, 3-particle frequencies are given by [See, for example, the review leading up to AI1 eq. (2.33b) ref. [6].] For QED, theq of a photon is zero, and so this formula becomes For the initial 3-particle evolution of the sequential diagrams of fig. 7(a-c) and virtual diagrams of fig. 8(h-k), this gives which is the frequency we quoted for xyyx in (4.12). For the sequential diagrams, the corresponding final 3-particle evolution has frequency for figs. 7(a,b) and its complex conjugate (Ω seq f ) * for fig. 7(c).

b. 4-particle evolution frequencies and modes
For 4-particle evolution, one just needs to repeat the derivation of AI1 section V.B [6], which studied the medium-averaged evolution of four high-energy gluons in (large-N c ) QCD. One can see from the diagrams of figs. 7-9 that, for large-N f QED, the only intermediate 4-particles states areēeēe, and so that is the only case we address here. (Beyond the large-N f limit, one would need to also considerēγeγ.) The only difference in the derivation is that the potential for four large-N c gluons [AI1 eq. (4.19) [6]] is replaced by the potential for (ē, e,ē, e), where b ij ≡ b i −b j and where the signs in front of the terms above are minus the product of the corresponding charges ±1. (E25) is algebraically equivalent to Proceeding as in AI1 [6], one finds the normal mode frequencies FIG. 40: QCD diagrams involving 4-gluon vertices, evaluated (for large N c ) in ACI4 [9]. In this figure, the solid lines above all represent gluons. As usual, the complex conjugates of the above interference diagrams should also be included by taking 2 Re[· · · ] of the above.
The corresponding eigenvectors C ± are which have been appropriately normalized so that We've chosen to work in the basis (C 41 , C 23 ) here, rather than the basis (C 34 , C 12 ) used in AI1 [6], in order to match the numbering used on sequential diagrams in ACI2 fig. 24 [7]. 35 One may check that a y a ⊤ y = (M ′ ) −1 , as implied by the normalization condition (E29).

Diagrams with instantaneous vertices
The calculation of the real double-spitting diagrams involving instantaneous vertices, shown in fig. 7(d-g), is very similar to the calculation of the QCD diagrams involving 4point gluon vertices, shown in fig. 40, which were computed in ACI4 [9].
a. The IĪ diagram Consider fig. 41(d), which is a labeled version of fig. 7(d). We will call this the IĪ diagram with "I" short for "instantaneous 4-fermion vertex." Its evaluation closely parallels that of 35 One may convert to the basis (C 34 , C 12 ) by simple permutation of the indices. That's because of charge conjugation symmetry, which means that the result for (ē, e,ē, e) is the same as that for (e,ē, e,ē), and then cyclically permute the indices of the latter to get back to (ē, e,ē, e) with (x 1 , x 2 , x 3 , x 4 ) → (x 2 , x 3 , x 4 , x 1 ).
the 44 gluon diagram of ACI4 section III.C [9], with the only change being that the overall factor coming from the two 4-gluon vertices in ACI4 eq. (3.12) [9] is replaced here by the factor To see this, compare the 4-gluon vertex rule of ACI4 fig. 10 [9] with the rule of fig. 39 here. From ACI4 eq. (3.16) [9] (times 3 to sum up the equal results from all three color routings), the QCD gluon result was The corresponding result here is then ∆t Ω + csc(Ω + ∆t).
Following our general procedure from AI1 [6] of subtracting out the vacuum pieces of each diagram (which must all cancel in the final result), this is ∆t x e y e z e (1 − x e ) 4 iΩ + ln 2. (E36) Adding this diagram to its complex conjugate gives what we labeled as the "(II)" contribution to the total answer (E1), x e y e z e (1 − x e ) 4 Re(iΩ + ) ln 2. (E37)

b. Diagrams with one instantaneous vertex
Now consider the interference diagram of fig. 41(e). Remember that our (large-N f ) convention is that x e is the momentum fraction of the particular final-state electron that is connected by an electron line to the initial electron.
(d) There are just a few differences with the similar calculation of the 4xȳ gluon process given in ACI4 section II [9]. First, comparing fig. 41(e) here with ACI4 fig. 12(a), our labeling conventions are a little different. The translation is that (x 1 ,x 2 ,x 3 ,x 4 ) = (−1, y, 1−x−y, x) in ref. [9] is permuted to (x 1 ,x 4 ,x 2 ,x 3 ) = (−1, x e , y e , 1−x e −y e ) here. Another is that the 4-point vertex factor of ACI4 eq. (2.9) for its color routing 4ȳx 2 should be replaced by Making these changes, and keeping track of the signs of the B ij in our fig. 37 rule, we find that ACI4 eq. (2.12) translates to The helicity sum analogous to ACI4 eqs. (2.14) is Using the QED formulas in appendix D here, we find that the initial-helicity average of (E41) above is ζ(x e , y e ) δnm (E42) Contrast to ACI4 eq. (2.16). Following ACI4, the final result, analogous to ACI4 eq. (2.25), is then where X seq x and Y seq x are as in (E13). Following the strategy of ACI4 section III.A, thexȳI diagram of fig. 41(f) here is like the mirror reflection of 41(e). Given how we have labeled particles, 37 this mirror reflection transformation can be achieved by the replacement and For the Ixȳ diagram, Ω seq f (x 1 ,x 2 ,x 3 ,x 4 ) is just a way to write the Ω seq f of (E24) in terms of the 4-particle variables (x 1 ,x 2 ,x 3 ,x 4 ). But the form (E46) has the advantage that, after the substitution (E45), it also gives the correct Ω for the initial 3-particle evolution forxȳI as depicted in fig. 41(f).
ThexIȳ diagram of fig. 7(g) vanishes by parity symmetry for the same reasons as thē y4x diagram in ACI4 section III.B.
We can summarize in a notation that parallels ACI4 [9] as closely as possible: 36 Similar to the discussion in footnote 33, absolute value signs have been judiciously included in (E43) so that ζ does the right thing under front-end transformations, which for the diagram at hand is implemented by (4.5). Note that the front-end transformation rule (4.5) for this diagram does not change the sign of x e y e z e , and so there is no need to write |x e y e z e | here like there was in (E16). 37 The labeling of figs. 41(e) and (f), and so the specifics of the transformation (E45), are slightly different than in ACI4 in order to maintain our convention here that x e always refers to the final-state electron whose electron line is connected to the initial-state electron.  From the starting formula (4.8), following the same steps as our earlier work [6,7] (AI1,ACI2) on real double splitting rates in (large-N c ) QCD, one obtains formulas that schematically have the same form. Namely, [which here must eventually be integrated over y e as in (4.7)]. In our QED application, Ω − = 0 (E27), and so above. The I new n are just like the I seq n of (E12), except that the (X, Y, Z) seq of (E13) are replaced by (X, Y, Z) new in which M seq f is replaced by M i (because the final 3-particle evolution in the xyyx diagram, fig. 29, involves the same particles as the initial stage of 3-particle evolution) and the a seq x are replaced by a y ′ = a y (because the particles that merge at the end of the interim 4-particle evolution are the same as the ones that split at its start): To study the UV divergence of xyyx, we want to identify which terms of the integrand in (F1) are important as ∆t → 0. For that, we need the small-∆t expansion of (F4). At leading order in ∆t, the expansion is the same for (X, Y, Z) new y and (X, Y, Z) new yy ′ : where the second equality uses (E30) and (E31). This approximation is fine for the Y 's and Z's in (F1), but it is inadequate for the combination X new y X new y ′ − (X new yy ′ ) 2 that appears in (F3), which is zero at the order of (F5).
Going to next order in ∆t, which gives One then finds that the I new n of (F3) are all O(∆t). Using these expansions in (F1) shows that all the terms in the integrand are finite as ∆t → 0 except for the term involving Z yy ′ I 1 , whose dependence on ∆t is This Z yy ′ I 1 term contains the small-∆t divergence represented by (4.25) in the main text.
[We will discuss the translation to (4.25) later.] Because (F7) blows up as (∆t) −2 , we also potentially have a sub-leading divergence (∆t) −1 if there are any corrections in the small-∆t expansion that are suppressed by only one more power of ∆t. Tracing back through (F7) and the expansions (F6), and noting that ln 1 − X 2 yy ′ /X y X y ′ = ln (X y X y ′ − X 2 yy ′ )/X y X y ′ in (F3b), the only correction that contributes at this order is the O (∆t) 0 correction to (F6e). Computing this correction requires the O(∆t) terms in (F6a) and (F6b).
There is a simplification that will allow us to work out a useful formula for the needed corrections to X new y and X new y ′ more generally than for the QED application of relevance here. So, even though Ω − = 0 for QED, let us be more general and rewrite (F4) as Noting that (4.9) and (E10) give x 1 x 4 (x 1 + x 4 )E = M i in (F5), the small-∆t expansions of the X's are then Then The expression a −1⊤ y Ω 2 a −1 y 11 sounds like a mess that depends on detailed formulas for Ω ± and a y . Happily, it can be recast into a very simple form.
b. Value of a −1 ⊤ Ω 2 a −1 As in appendix E 2 b above, return again to the derivation of eigenfrequencies Ω ± and normal modes in AI1 section V.B [6]. From AI1 (5.17), the relevant Lagrangian has the form where C is (C 41 , C 23 ) or whatever permutation you want, and M and K are 2 × 2 matrices in that basis. The 1 2 C ⊤ K C above encodes the potential V of the double harmonic oscillator problem, e.g. (E26) in our case. For our basis choice (C 41 , C 23 ), M is given by (E30).
The equation of motion from (F12) is and the corresponding normal-mode eigenvalue problem is We can rewrite both cases (±) of this equation simultaneously as where is the 2 × 2 matrix a y [(E14) in this paper, or whatever permutation is relevant to one's choice of basis for C]. Multiplying (F15) by a ⊤ on the left, Because the normal modes are orthogonal with respect to M and then normalized (E29) so that we have a ⊤ Ma = 1 and so This result is independent of details about the eigenfrequencies or eigenmodes and which basis we pick for C.
In the application to xyyx, we want in particular the first element of this matrix, in the (C 41 , C 23 ) basis (numbered as in the right-hand diagram of fig. 29). In that context, K 11 corresponds to the (complex) spring constant of the problem if we were to set C 23 to zero. That would be the same as setting b 2 = b 3 and so the same as placing the e andē of the photon self-energy loop on top of each other. In that case, the eē pair is, with regard to charge, indistinguishable from the photon that created them (i.e. no charge in the QED case here). In this case, the potential in the 4-particle evolution of fig. 29 would be the same as that in the preceding 3-particle portion. But that means that K 11 must be identical to the spring constant for the initial 3-particle evolution, and so in the context of the xyyx diagram. One may verify this general relation in our QED case by starting from the potential (E26), following the method of AI1 section V.B [6] to construct in the (C 41 , C 23 ) basis, and then checking (F22) using (4.9), (4.12), and (4.22). A similar check works in the case of (large-N c ) QCD.

c. Putting it together
Using the results above, (F11) becomes Using this and our earlier expansions in (F1) gives the unregulated (ǫ=0) divergence This is the origin of (4.11) in the main text. As we will see below, theĪ unregulated above is the unregulated (ǫ=0) version of the I introduced in (4.30) of the main text but generalized here to be appropriate for any sign of M i . 38 We will see that the X of (4.25) turn out to represent the X new of (F4) without the |M i |Ω i terms, i.e.
The expansions (4.26) of the X then follow from the expansions (F10) of the X. Finally, we note that the relationship (E19) and the explicit values (4.9) of the x i can be used to rewrite (F25) as The overall factors of P e→e P γ→e are the source of the similar overall factors in (4.25), except that we will need to generalize the derivation to d = 2 − ǫ transverse dimensions.
38 See footnote 36. The absolute value signs on all |M i | appearing in (4.25), (F10) and later (F31) are accounted for by |M i |Ω i = M i Ω i sgn(M i ) = M iΩi . The notationΩ defined by (F28) is unrelated to the notationΩ f in AI1 section VI.B [6]: There are only so many accent marks, and we've been forced to recycle.

Needed generalizations to d=2−ǫ
In the case of QCD g → ggg, the dimensional regularization of diagrams was carried out in ACI3 [8]. We can adapt intermediate results from that paper if we (i) first look at its sequential diagram result for xyxȳ [ACI3 eq. (5.10)], (ii) convert to QED by making the same modifications as in section F 1 of this paper, and then (iii) adapt the results to the virtual xyyx diagram as in (4.8). But, for the reasons described in section IV C 3, we also need to backtrack in the derivation of ACI3 and never expand the Bessel functions K d/4 ( 1 2 |M|ΩB 2 ) [ACI3 eq. (4. 15)]. The result is the following generalization of the divergent Z yy ′ I 1 term of (F1) to d dimensions: where the M i , Ω i and (X , Y, Z) are defined exactly the same as for d=2.
To make contact with the derivation in ACI3, compare (F31) above to the Z term in ACI3 (4.14) for the QCD xyȳx diagram. The structure is the same. In addition to a factor of 2 related to converting the QCD group factors for xyȳx to QED, and the different way (ᾱ,β,γ) are contracted here to make dᾱ+β+γ for the Z term of xyyx, there is one minor change in the overall pre-factor concerning the powers of µ and E in (F31) above. µ is the renormalization scale, which we introduce by writing the d-dimensional coupling constant g d , which has dimensions of (mass) ǫ/2 , as where g is the dimensionless coupling constant for ǫ=0. As a result, the dimensionless α 2 EM in (F31) is associated with a factor of µ 2ǫ , as written explicitly in the pre-factor. The overall power of E then follows from dimensional analysis if one takes the convention, as in ACI3 [8], that the d-dimensional generalizations of (ᾱ,β,γ) are defined to be dimensionless. 39 39 To wit, there was a mistake in the prefactors in ACI3 [8], which are off by an overall factor of (µ/E) 2ǫ .
This mistake did not matter there because the 1/ǫ poles all cancel when one adds up all of the real-double splitting diagrams (because there should be no UV divergence in the total result if there are no loops in the amplitude or in the conjugate amplitude). In that case, the ǫ dependence of an overall factor common to all diagrams will not matter when we set ǫ=0 at the end of the day. In this paper, however, the 1/ǫ poles do not cancel when we sum the virtual diagrams of fig. 8-they can't, because we know we need to get coupling renormalization at this order. ACI3's error in the overall factor comes in the paragraph of ACI3 appendix A concerning ACI3 eq. (4.3), which forgets that the d-dimensional coupling is dimensionful.
We may now use the d-dimensional generalization [8] 40 e→e (x e ) P We will see later that we do not need the d-dimensional version of the DGLAP splitting function P (d) e→e (x e ) in (4.25) because it is common to xyyx and the leading-order process xx. But we will need the d-dimensional version of the other DGLAP splitting function P which reproduces the usual result in the case ǫ=0.

The subtraction D 2 (∆t)
We now turn to the subtraction D 2 (∆t) introduced in (4.16) and (4.17), which will allow us to (i) turn the d=2 expression (F1) for xyyx into a convergent integral that can be done numerically, corresponding to the d(∆t) [F 2 (∆t) − D 2 (∆t)] term above, and (ii) cancel the a dependence of the ∆t<a contribution (F34), as in the first two terms above. We will choose D 2 (∆t) proportional to (4.18). To get the ultimate proportionality 40 The argument for this relation is the same as that for ACI3 (5.17) [8]. However, fixing the overall normalization error described in footnote 39 of the current paper modifies the (1 − x) d−1 in the second denominator of that equation to 1 − x. Like the overall normalization, this correction does not affect any of ACI3's final results because of the cancellation of divergence there, but it is important for our current work. 41 See, for example, eq. (16) of ref. [32], which one may verify independently. Our ǫ is their 2ǫ, and their T R is 1 in the QED case we consider here. The result implicitly depends on the convention [33] that the trace of the Dirac identity matrix is simply defined as tr(1 Dirac ) ≡ 4 in d dimensions, which is part of Conventional Dimensional Regularization for fermions.
constant, see (F25-F27), but first we will define a D Combining this with the result (4.31) for I, we see that the a dependence cancels, leavinḡ (F39) Now multiply this by the prefactors shown in (F34) to convert I into the ∆t<a result for xyyx to get 42 e→e (x e ) P We'll leave it in this form for the moment.
As designed, this is an integral that can be done numerically.

Integration over y e
For the virtual diagram, we need to integrate over y e as in (4.7). For the subtracted piece (F41), this is another integral we will do numerically: (F42) For the other piece (F40), we will do the y e integral analytically. We start by changing integration variable from y e to y e ≡ y e /(1 − x e ): To carry out renormalization as in (4.37), we need the leading-order result [dΓ/dx e ] xx in d transverse dimensions. That can be taken from the similar result in ACI3 [8] for g → gg except that C A α s P g→gg (x) in QCD is replaced by α EM P e→e (x e ) here: 43 where B(x, y) ≡ Γ(x) Γ(y)/Γ(x+y) is the Euler Beta function and the factor of µ ǫ associated with α EM comes from (F32 is given by (F41). Eq. (F47) is the result (4.38) quoted in the main text, but generalized here to handle either sign of M i and so handle front-end transformations such as in (4.39). Since the renormalized result (4.38) is finite, we may use the d=2 version (B5) of (F45) there, which is So, as previously promised, we never need the d-dimensional version of the structure function factor P (d) e→e (x e ) that was common to the leading-order xx and the virtual correction xyyx.

7.q and dimensional regularization
Throughout this paper, we have usedq as an independent parameter describing interactions with the medium, which could either be calculated theoretically (in certain limiting cases) or used as a phenomenological parameter. One might worry whether or not the use of dimensional regularization requires knowing the O(ǫ) corrections to the 3+1 dimensional value ofq, given that various of our formulas along the way have involved 1/ǫ divergences multiplying expressions that depend onq through complex frequencies Ω. Fortunately, this is not an issue. Imagine that we used the full d-dimensional value ofq (whatever it is) everywhere in our intermediate calculations. Since our final, renormalized results are finite expressions where the ǫ → 0 limit is taken, theq in those expressions can then be replaced in the last step by 3+1 dimensionalq, just as we did for P e→e (x e ) above. 44 where the X have the small-∆t expansions (4.26) X y = X y ′ = − iM ∆t + iMΩ 2 ∆t 3 + · · · , (G2a) Plugging these expansions into the exponential of (G1) gives where the last exponential factor represents the corrections from the second terms in the expansions (G2) and O(· · · ) means "of order." Because of the bound (4.29) on the sizes of both B and B ′ that contribute to I, the exponent of this second exponential is small: We will need this correction, but we can treat it as perturbative, rewriting (G3) as (G5) Formally, it will be convenient to implement this expansion by introducing a redundant parameter ξ=1 into (G2), and then think of I as a function I(ξ) of ξ. The above discussion then translates to I = I(ξ) ξ=1 = I(0) + ξ I ′ (0) + ξ 2 2! I ′′ (0) + · · · ξ=1 = I(0) + I ′ (0) + O(a).
For the same reason that the unregulated (ǫ=0) discussion of appendix F 2 did not require the yet-higher order terms not explicitly shown in (G6), we will not need them here either. They give vanishing contribution to I for a → 0.

Units
In order to simplify the presentation of calculations, in this appendix we will evaluate I in units where M = 1 and Ω = 1.
At the end, we will be able to put M and the complex-valued Ω back into the answer using the scaling properties of the definition (G1) of I and then analytic continuation of Ω back to complex values.
So (G1) for I becomes

Doing the B integrals
To do the B integrals above, it is convenient to take advantage of the form of the exponential to replace The B integral of the exponential gives B,B ′ exp − 1 2 (X y + p)B 2 − 1 2 (X y ′ + q)B ′ 2 + X yy ′ B · B ′ = (2π) d det X y + p −X yy ′ −X yy ′ X y ′ + q −d/2 = (2π) d (X y X y ′ − X 2 yy ′ ) + X y ′ p + X y q + pq Now use (G11), giving . (G13) Now specialize to our case by using (G6) to get As discussed earlier, the relative O (∆t) 2 corrections can be dropped because they will lead to a convergent ∆t integral in d=2 that vanishes as a → 0. Using (G14) in (G10) gives ∞ 1 dp dq (p 2 − 1) −ǫ/4 (q 2 − 1) −ǫ/4 × a 0 d(∆t) (∆t) (d+2)/2 −i(p + q) + (pq + ξ)∆t −(d+2)/2 + O(a). (G15) We find it convenient to change integration variables to to get 46 where F = 2 F 1 is the hypergeometric function. In our application (G17), a 0 d(∆t) (∆t) (d+2)/2 −i(u + v) + (1 + ξuv)∆t −(d+2)/2 (not to be confused with any other use of the letters α or β in this paper, but there are only so many letters in the alphabet). Use the hypergeometric transformation together with F (a, 0; c; z) = 1 [which follows from the series expansion that defines the hypergeometric function], to rewrite (G19) as a 0 d(∆t) (∆t) (d+2)/2 −i(u + v) + (1 + ξuv)∆t (G22) The advantage of this rewriting is that the second term on the right-hand side of (G22) is finite if we set d = 2. That is, the 1/ǫ divergence arising from the ∆t integration is isolated in the first term. Moreover, we will see that the integration of the finite second term over (u, v) is also finite, so we do not need to keep dimensional regularization in order to do those (u, v) integrals: we can just set d=2 there and be done with it. All together, using (G22) in the expression (G17) for I gives with (dimensionally regularized) divergent piece The ignorability of higher-order terms in the ξ-expansion (G27) of I in the a→0 limit will hold for the total I = I div + I regular but turns out not to hold separately for I div and I regular . A check we should make on our calculation is that the O(ξ 2 , ξ 3 , · · · ) terms indeed cancel between I div and I regular . For this purpose, it will be sufficient to check that I ′′ div (ξ) + I ′′ regular (ξ) = O(a) for arbitrary ξ ≤ 1. So, for later reference, we note here that differentiating the integrand of (G25) twice with respect to ξ gives 7. Evaluating I div Rewrite (G24) as As before, we take the ξ expansion (G27) of I and so here of A. The integral does not converge for d = 2−ǫ near 2, and so we cannot simply expand the integrand in powers of ǫ. Dimensional regularization tells us to imagine doing the integral in dimensions d < 1 where it converges and then analytically continuing the result in d. But we can make that job easier if we first rewrite (G34) as