Strong- vs. weak-coupling pictures of jet quenching: a dry run using QED

High-energy partons ($E \gg T$) traveling through a quark-gluon plasma lose energy by splitting via bremsstrahlung and pair production. Regardless of whether or not the quark-gluon plasma itself is strongly coupled, an important question lying at the heart of philosophically different approaches to energy loss is whether the high-energy partons of an in-medium shower can be thought of as a collection of individual particles, or whether their coupling to each other is also so strong that a description as high-energy `particles' is inappropriate. We discuss some possible theorists' tests of this question for simple situations (e.g. an infinite, non-expanding plasma) using thought experiments and first-principles quantum field theory calculations (with some simplifying approximations). The physics of in-medium showers is substantially affected by the Landau-Pomeranchuk-Midgal (LPM) effect, and our proposed tests require use of what might be called `next-to-leading order' LPM results, which account for quantum interference between consecutive splittings. The complete set of such results is not yet available for QCD but is already available for the theory of large-$N_f$ QED. We therefore use large-$N_f$ QED as an example, presenting numerical results as a function of $N_f\alpha$, where $\alpha$ is the strength of the coupling at the relevant high-energy scale characterizing splittings of the high-energy particles.


A. Overview
Consider the cartoon in fig. 1a depicting energy loss of a high-energy parton traversing a lower-energy quark-gluon plasma (QGP) produced in a relativistic heavy-ion collision. This figure and this paper focuses on the in-medium evolution, ignoring the complications of both (i) early-time initial vacuum-like radiation associated with the hard parton-level collision that scattered the high-energy parton out of the beam direction in the first place and (ii) the late-time hadronization when a high-energy parton leaves the quark-gluon plasma. Generically, energy loss of high-energy relativistic particles traveling through a medium is dominated by bremsstrahlung and pair production, as depicted by the splittings in the cartoon. In detail, this picture of energy loss implicitly assumes that the high-energy particle and its high-energy daughters and granddaughters et cetera may be thought of as individual particles, which in quantum field theory implies a perturbative description of the states of the high-energy particles (even if the underlying quark-gluon plasma is strongly coupled). Such a perturbative treatment requires that the running QCD coupling α s (µ) be "small" at the relevant scale µ that characterizes the high-energy splittings. This is in contrast to the separate question of whether the quark-gluon plasma itself is weakly-or strongly-coupled, which depends on the strength of α s at QGP scales, such as α s (T ) or α s (ξ −1 ), where T is the temperature and ξ is the QGP chromo-electric screening length. (Initial vacuum-like radiation associated with the hard event that created the initial particle is not shown. Vacuum-like radiation and then hadronization of a high-energy particle that leaves the medium is also not shown.) (b) Depiction of a branch that deposits all of its energy in the medium.
Though real-world quark-gluon plasmas are generally considered to be strongly coupled (or to at the very least require complicated reorganizations of small-coupling expansions [1]), the situation has been less clear with regard to the coupling of each high-energy parton to its daughters in the energy-loss picture of fig. 1. On the "weakly-coupled splitting" side, there is a large literature built upon the formalism of Baier et al. [2,3] and Zakharov [4] (BDMPS-Z), which treats high-energy partons as individual particles. At high energy, the dependence of splitting on medium properties may (with caveats) for many purposes be summarized by the medium parameterq, whether the medium itself is weakly or strongly coupled.q is the proportionality constant in the relationship Q 2 ⊥ =q ∆z, where Q ⊥ is the transverse momentum that a high-energy parton picks up, relative to its initial direction of motion, after passing through length ∆z of the medium. Alternatively, on the "stronglycoupled splitting" side, where high-energy partons cannot be treated individually, there are studies of how energy loss behaves at extremely large coupling in QCD-like theories with gravity duals [5][6][7]. 1 Because of all the complexities that enter describing the many different stages of relativistic heavy-ion collisions, it has been difficult to settle the issue of weakly vs. strongly coupled splitting from experimental data. So it would be useful to have theoretical calculations we could do, for some simplified thought experiments, that would shed light on whether we are in the weakly-coupled or strongly-coupled splitting limit for energy loss.
As we shall review, previous authors [9][10][11] have made leading-log studies of the correction to BDMPS-Z splitting rates due to additional, softer gluon bremsstrahlung happening simultaneously with the underlying splitting process. They found that such corrections are large but are also universal in the sense that they are completely absorbed (at leading log level) by accounting for similar corrections toq. One may then pose a refined question about the merits of a weakly-coupled splitting picture of in-medium showering: how large are the corrections to this picture that cannot be absorbed into an effective value ofq? 2 In this paper, we will give some examples of measures that can be used in theory thought experiments to address this question, and we will discuss advantages and disadvantages of those measures. Though our ultimate motivation is to eventually address the question for QCD, the full set of QCD calculations required is not yet available. Here we will show concrete results for a different theory: large-N f QED. For that case, we give quantitative measures of the size of corrections to the weakly-coupled splitting picture as a function of the size of N f α EM .
For the sake of simplicity, we are going to focus on high-energy particles which completely stop and deposit their energy in the medium, like the first branch of fig. 1a, isolated in fig.  1b. Phenomenologically, one could view this as a choice to focus on partons whose energy is large compared to that of typical plasma particles but not so large that they punch all the way through the plasma. Alternatively, one could instead consider arbitrarily high energy partons but imagine, as a theory exercise, that the extent of the plasma were large enough for them to stop completely. Whether in-medium showers of high-energy particles have weakly-coupled or strongly-coupled splittings in this case is an interesting and important 1 For an example of a phenomenological model that attempts to combine certain aspects of both the weakly-coupled and strongly-coupled pictures of splitting, see the hybrid model of ref. [8]. 2 There is precedence for a similar distinction in a very different context: For calculations of a weakly-selfinteracting quark gluon plasma [that is, α s (gT ) small, which we do not assume here], calculations of the shear viscosity to entropy ratio η/s of the quark-gluon plasma find that the next-to-leading order (NLO) correction in α s (gT ) is very large but almost completely accounted for by the also-large NLO correction toq [12]. (Technically, the calculation covers only "almost" all next-to-leading order contributions to η/s [12]. ) We should emphasize that their expansion in the quark-gluon plasma coupling is not our expansion in the high-energy splitting coupling, and what they call a NLO correction toq is something that in our expansions will already be part of the leading-order value ofq, which we will later refer to asq (0) . first question. And studying that limit provides a way to examine the issue independent of the complicated details necessary to discuss all cases relevant to the full phenomenology of energy loss.
Throughout, we will use the phrase "quark-gluon plasma" as a generic term for the QCD background generically produced in heavy-ion collisions without committing to any more specific, detailed picture of that background other than that scattering of high-energy particles from the medium can be characterized by some value ofq.

B. Overlapping Formation Times
In any single example of a shower, the exact time or place where each splitting occurred has some ambiguity, and the extent of that uncertainty is known as the formation time or formation length for that splitting. In terms of calculations of the splitting rate, consider multiplying an amplitude times a conjugate amplitude to get a rate. The formation time corresponds, parametrically, to the largest time difference |t −t| for which a splitting at time t in the amplitude non-negligibly interferes with the same splitting at timet in the conjugate amplitude, as depicted in fig. 2. At high energy, the corresponding formation length in the direction of motion is the same as the formation time (times the speed of light), and they grow with energy. For example, for bremsstrahlung with energy xE from a parton with energy E, the formation length is parametrically 3 in the limit of a thick medium. When the formation time becomes large compared to the mean free path for collisions with the medium, as in fig. 2, it causes a reduction in the bremsstrahlung rate known as the Landau-Pomeranchuk-Migdal (LPM) effect [14,15]. A similar reduction occurs for pair production. In fig. 3a, we show a shower where we have qualitatively depicted the formation lengths for each splitting by ovals. This picture implicitly assumes a hierarchy of scales in which the distance between consecutive splittings is large compared to the formation lengths for each splitting, so that the probability of each splitting is independent. If true, one may take as a starting point for simulating this type of in-medium shower a calculation (or model) of single-splitting rates. As simulation time progresses, one could roll dice every time increment to determine whether a new splitting happened in that time interval. Such a model of independent dice rolls for each splitting assumes that the chance of two or more consecutive splittings having overlapping formation times, such as depicted in fig. 3b, is small.
Consider for a moment the case of nearly-democratic splittings, where the two daughters of each splitting have roughly comparable energy. In general, at very high energy, each 3 For a very brief review of formation times in the language of this paper, see section II.A.1 of ref. [13]. The QED formation time essentially goes back to Migdal [15], though he did not use this language. (qt form is proportional to what Migdal would call m 2 /s for bremsstrahlung and m 2 /s for pair production.) Interference between high-energy bremsstrahlung at time t in the amplitude and timē t in the conjugate amplitude, as the high-energy particles repeatedly scatter (represented by the brown dots) by small angles from the medium. In QCD, the particle represented by the non-curly line could be a quark or a gluon. For photon bremsstrahlung, photon scattering from the medium can be ignored. formation length of scattering from the medium offers one chance for splitting, and the chance of such a splitting is parametrically of order the coupling α associated with the splitting vertex. The typical distance l rad between consecutive splittings is then If the relevant value of α is small enough, this leads to the hierarchy of scales shown in fig.  4, leading to the picture of independent splittings in fig. 3a. The probability of any two particular splittings overlapping, like the two in fig. 3b, would be O(α).
Since α runs and depends on scale, we need to identify the relevant renormalization scale for the α in (1.2). That α characterizes the strength of the vertex where the high-energy parent splits into its high-energy daughters, and the scale we need is therefore the distance scale characterizing the separation of the daughters during the splitting process, i.e. during the formation time. That separation b is transverse and so is related by the uncertainty principle Q ⊥ ∼ 1/b to the relative transverse momentum Q ⊥ ∼ √q t form that the high-energy particles pick up during the splitting process. So we are interested in α(Q ⊥ ), with (1.1) giving Q ⊥ ∼ (qE) 1/4 . The above argument that (1.3) is the relevant scale for α is qualitative, but it has been checked by explicit calculation in the case of large-N f QED [13]. Note that though this scale grows with energy E, it does not grow quickly! For example, forq ∼ 1 GeV 3 for quark-gluon plasmas and E as large as ∼ 100 GeV, the parametric estimate (1.3) is Q ⊥ ∼ only a few GeV (in the theoretical limit of a thick medium). If α(Q ⊥ ) is too big, then splittings will usually overlap, and we cannot then use a weaklycoupled description of splittings. A way to diagnose this problem is to first formally assume α is small, calculate the corrections to shower development due to overlapping formation times, and then see how large those corrections are when one evaluates them for realistic values of α. If the corrections are of order 20%, for example, then the basic picture of weaksplitting behind figs. 1 and 3a is reasonable. If the corrections are 200%, then it is not, and the number of high-energy particles present at any moment in the shower would not (even approximately) be a well-defined concept.
Throughout this paper, the coupling α will generally refer to α(Q ⊥ ) [and not to α(T )].

C. Stopping distances
Hard bremsstrahlung is more efficient at decreasing the energy of particles in the shower than soft bremsstrahlung. In more general language, nearly-democratic splittings are more efficient for development of the shower than splittings where one of the daughters is soft. To study the size of overlap effects on the development and stopping of showers, we would like to find simple quantities that naturally weight different processes by how much they contribute to the degradation of the energies of the particles in the shower. A natural quantity to consider is simply the characteristic length of the shower to where it ends and deposits its energy into the medium. We'll generically refer to this as the average "stopping distance" ℓ stop .
There are different variants of stopping length one could define. If the original high-energy particle carries a conserved charge, like fermion number, then one could track the statistical distribution ρ(z) of the charge density deposited in the medium as a function of distance z from the initial position of the initial high-energy particle. 4 (Transverse displacements will be small compared to z, so we need not distinguish between distance and longitudinal distance.) Following ref. [16], we then define the average "charge stopping distance" as the 4 Once deposited, charge will then flow with the medium and also diffuse hydrodynamically. Our definition here of ρ(z) does not follow that evolution but instead refers to where the charge was deposited when it became part of the medium. first moment of the distribution: where q 0 is the charge of the initial particle (e.g. q 0 =1). Alternatively, we could do the same with the energy density ε(z) deposited in the medium to define an "energy stopping distance," In principle, the very last stage for particles from the shower stopping in the mediumwhen their energies finally become comparable to the characteristic energy scale T of the plasma-is not described by the high-energy approximations that we shall use throughout this paper. Fortunately, the effect of that last stage has only a parametrically small effect on the stopping length for showers initiated by high-energy particles E ≫ T [16]. To see this, consider the scale (1.2) characterizing the distance between consecutive splittings. For nearly-democratic splittings, this is As a crude estimate, imagine that the energies of individual particles drop by a factor of 2 with each splitting. Then the total distance the shower would propagate would be This is a convergent series giving, parametrically, The error in the treatment of stopping when the particles degrade to energy of order T is just a piece of order l rad (T ) ∼ α −1 T /q in the series (1.6), which is suppressed by T /E compared to the sum (1.7). We could now ask what effect overlapping formation times have on these stopping lengths. Is it a small or large effect for relevant sizes of the coupling α? Consider an expansion (1.8) Here ℓ stop is the result in the approximation that no splittings interfere with each other, in which case the statistical development of the shower can be computed based just on single-splitting probabilities. ∆ℓ stop is the O(αℓ stop ) correction due to overlapping formation lengths, computed in the formal limit that α is small. We might then consider taking the ratio ∆ℓ stop as our measure of whether splitting in showers is weakly or strongly coupled for a given size of α.
There is a problem with this measure for QCD. In the parametric estimates leading to (1.9), we have so far only discussed nearly-democratic splitting. However, soft gluon bremsstrahlung (x≪1) is associated with shorter formation lengths than hard gluon bremsstrahlung, as reflected in the small x dependence of (1.1a). That means that soft gluon bremsstrahlung is less LPM suppressed and so happens more often than hard gluon bremsstrahlung. Various authors [9][10][11] have found that overlap effects of soft-gluon bremsstrahlung during harder-gluon bremsstrahlung, as depicted in fig. 5, correct the harder emission probability by an amount that is formally suppressed by order 5 α s ln 2 L τ 0 (1.10) instead of α s , where τ 0 characterizes the scale of the mean-free path for scattering from the medium and reflects a characteristic scale of the medium. (Here we write α s instead of α because this part of the discussion is specific to QCD and does not apply to QED.) In the context of a thick medium, as considered in this paper, the L in the logarithm should be interpreted as the formation time for the harder splitting. If the harder splitting is nearly democratic, then L ∼ E/q and so the effect of a soft emission during a nearly-democratic splitting is suppressed by order which is no suppression at all for energy large enough that the double log compensates for the smallness of α s (Q ⊥ ). This effect on splitting rates will result in the same soft double-log enhancement to the ratio considered in (1.9), so that (1.12) However, the double-log corrections to splitting rates have a very special form [9][10][11]: they can be absorbed into a redefinition of the effective value ofq, exactly corresponding to the result of an earlier calculation [18] of the double-log soft-bremsstrahlung correction δq to the definition ofq itself. At first order in α, in leading-log approximation, they found where R indicates the color representation of the particle whose deflection is being described byq, where C A =N c is the quadratic Casimir for the adjoint color representation (i.e. for a bremsstrahlung gluon), and where ≈ is our notation for leading-log approximation. Moreover, the leading logs can be summed to all orders in α s [18] with the result that the effective value ofq scales with L asq eff ∼ L γ [9] for sufficiently large values of L, with γ = 2 C A α s /π. For splitting in the limit of a thick medium, this corresponds to the scalinĝ (1.14) for the effectiveq to use in calculations of nearly-democratic splitting. If inserted into the parametric formula (1.7) for the stopping length, this gives 15) up to yet higher order corrections in the exponent. This is an important result for theory because it explains how, in QCD-like theories where strong-coupling results are also known (namely theories that have gravity duals), the weak-coupling behavior ℓ stop ∝ E 1/2 can change as α s is increased to move toward the known strong-coupling behavior ℓ stop ∝ E 1/3 [5][6][7]. The above E −γ/4 relative correction from soft overlapping emissions to the weak-splitting result E 1/2 is a parametrically large effect at large E, even if γ were small. Moreover, the relative correction γ/2 to the exponent is controlled by the size of √ α s rather than (the parametrically smaller) α s . The upshot is that it has already been known for some time that corrections due to overlapping formation lengths are very significant in the case of overlap with soft emissions, but that this effect can be absorbed into an effective value ofq. On the more phenomenological side, Zakharov [20] has recently calculated that other, non-logarithmic corrections toq at relative order α s are comparably important in applications of interest and can change the sign of the effect.
In this paper, we want to address the question of how large are the corrections due to overlapping formation times that cannot simply be absorbed into an effective value ofq? There are two ways to do this. One is to look for some different measure that is insensitive to the size ofq in the first place and so avoids double-log enhanced corrections. We'll propose a (partly) successful choice in section I D below.
Alternatively, one can re-organize the expansion (1.8) so that theq used in the leading term is calculated using the effective valueq eff ofq instead of the original valueq (0) , wherê q eff =q (0) + δq when formally expanded to first order in α. There is an ambiguity in exactly how one definesq eff in this context, however, because the identification of the scale L in (1.13) as the formation length ∼ E/q is only a parametric identification of scale. Similar to ambiguities associated with choosing renormalization scale in perturbation theory, or factorization scales for defining parton distribution functions at next-to-leading order, the exact choice of L here is a matter of convention as long as one chooses L ∼ E/q in order to eliminate the large double logarithm in the re-organized small-α expansion. (There are further concerns one can have about sub-leading, single log terms, but we leave that for later discussion.) Though we have motivated our discussion with QCD applications, in this paper we are going to give example calculations for the case of large-N f QED, where N f is the number of electron flavors. For this theory, all the overlap corrections to BDMPS-Z splitting rates that will be needed are already available [13]. (The only reason for the large-N f limit was to reduce the number of diagrams that needed to be calculated in ref. [13].) The behavior of the LPM effect in QED is qualitatively similar to the LPM effect in QCD for the case of nearlydemocratic bremsstrahlung or pair production but is qualitatively different for the case of softer bremsstrahlung. In QCD there is less LPM suppression of softer bremsstrahlung; in QED there is more. In particular, QED does not have the double-log enhancement (1.12). It does, however, have a different type of logarithm (which QCD also has): for any quantity that depends on coupling α at leading order, the next-to-leading-order (NLO) correction will necessarily have logarithmic dependence on the choice of renormalization scale. As we'll discuss, this logarithm has the form (1.16) As discussed qualitatively before in the context of (1.3), the renormalization scale should be chosen of order Q ⊥ ∼ (qE) 1/4 , here in order to eliminate large logarithms (1.16) in the expansion of ℓ stop in α. But, like our discussion of L for absorbing QCD double logs, the exact choice to make for µ is ambiguous. Since our goal is to judge the relative size of NLO corrections, we must therefore account for this ambiguity by considering a reasonable range of µ.
We will comment at the end of the paper on what one might take away from this and other QED results.

D. Shape of stopping distribution
Ideally, it would be nice if in the QCD case we could sidestep the scale ambiguities associated with absorbing the corrections toq and instead construct a thought-experiment observable that is independent of the size ofq. We have a (partly) successful proposal (we'll explain the caveat "partly" later on) which has the added benefit of pushing the renormalization-scale ambiguity of α(µ) off to next-to-next-to-leading order (NNLO).
The average stopping distance ℓ stop was the first moment (1.4) of the spatial distribution ρ(z) or ε(z) of where charge or energy is deposited by the shower. But one could also ask after the width σ of the region across which those quantities are deposited, as depicted in fig. 7. As we'll see later, both σ and ℓ stop are of the same order: So the scale ofq will cancel if we consider their ratio σ/ℓ stop . We propose computing the dimensionless ratio, which could be taken for either the deposited charge or deposited energy distributions, similar to (1.4). We will discuss how to compute the formal expansion of (1.19) in powers of α(Q ⊥ ), writing (1.20) Here, as in (1.8), the superscript "(0)" indicates an approximation based only on the results for single splitting rates, ignoring overlapping formation length effects. χα represents the relative size of the effect of overlapping formation lengths on σ/ℓ stop , calculated to first order in α. Our proposed test is whether or not χα is at least somewhat small compared to 1. One may also look at dimensionless ratios involving higher moments of the distributions, but for the moment we focus on (1.19). Note from (1.18) that the leading-order dependence on α will also cancel in the ratio σ/ℓ stop . Since the leading-order term (σ/ℓ stop ) (0) in (1.20) is independent of α = α(µ), there will be no explicit ln µ dependence in the NLO correction. So there will be no explicit renormalization scale ambiguity when we discuss the relative size χα of the NLO correction as a function of α. In contrast, we had to deal with the ambiguity when discussing the relative size of the NLO correction in the expansion (1.8) of ℓ stop because there the leadingorder result ℓ (0) stop depended on α. As a concrete example, a preview of our results in large-N f QED of the relative size χα of the NLO correction to σ/ℓ stop is The stopping distance and width characteristic of where charge or energy is deposited by a shower. (The above curve is presented here for qualitative purposes, but the shape depicted just happens to be the precise leading-order numerical result for QED charge deposition ρ(z)see appendix B. In that case, the specific "∼ σ" line shown above has length 2σ.) The relative size χα of the NLO correction in large-N f QED to the ratio σ/ℓ stop for the deposited charge distribution. The horizontal axis is N f α, where it is not necessary for this ratio at this order to distinguish the precise choice of renormalization scale (or scheme) other than that µ ∼ (qE) 1/4 parametrically.
We plot this in fig. 8 for the sake of visual comparison to fig. 6. The two measures give roughly similar conclusions for large-N f QED: For N f α ≃ 1, overlap corrections are large, in which case the weakly-coupled splitting picture (as in fig. 1) of the high-energy particles in showers would not be a very good picture. In contrast, for N f α ≤ 0.2, the overlap corrections are only a modest effect. It was of course a foregone conclusion that parametrically the weakly-coupled splitting picture would break down at N f α ∼ 1, but it was not previously clear whether quantitatively that meant N f α ≃ 1 or N f α/2π ≃ 1 or something else.

E. Why not talk about dE/dz?
It is more common and traditional to package discussions of energy loss into the evaluation of the rate dE/dz of energy loss of a particular particle per unit length of medium. So why have we instead been focusing on stopping lengths and the charge deposition distribution? As we'll explain in section II B, the definition of dE/dz becomes ambiguous, in the general case, once one considers corrections from overlapping formation times. In this paper, our philosophy is to consider tests of the importance of overlapping formation times which can be applied in general situations.
That said, it will turn out that dE/dz is unambiguously defined for the example we use for numerics in this paper: electron energy loss in large-N f QED. (It is also well defined for quark energy loss in large-N c QCD.) So for this case, we can make contact with the more traditional dE/dz by showing in figure 9 a plot of the relative size ∆(dE/dz)/(dE/dz) (0) of the NLO correction to dE/dz, analogous to our plot in fig. 6 for ∆ℓ stop /ℓ (0) stop .

F. Outline
In the next section, we first state a little more clearly some of our assumptions. Then we present in the main text the simplest derivations (that will be adequate for the explicit case of large-N f QED analyzed in this paper) of the moments z n of the charge deposition distribution ρ(z). From those results, we derive explicit integral formulas for ∆ℓ stop /ℓ stop and the measure χα in terms of leading-order and NLO contributions to differential rates dΓ/dx for splitting. We leave discussion of energy stopping distances, and more general techniques that can handle QCD double logarithms, to appendices. In section III, we then discuss in more detail the origin of the leading-order and NLO contributions to the relevant rates dΓ/dx, making contact with the explicit results and expressions given for large-N f QED rates derived in ref. [13], which are then used to produce our final numerical results of figs. 6 and 8. The NLO contributions to the splitting rates contain both (i) overlap corrections ∆Γ/dx dy for two consecutive emissions (like the overlap shown in fig. 3b), and (ii) corresponding NLO virtual corrections to leading-order, single-splitting results. At the moment, (ii) has not yet been computed for (large-N c ) QCD but has been for large-N f QED. In section IV, we discuss more about renormalization scale dependence and discuss how the breakdown of QED at high momentum scales (the Landau pole) is not a concern even for the relatively large values of N f α EM we have investigated.
In large part, we motivated consideration of the relative size χα of the NLO correction to σ/ℓ stop ( fig. 8) by looking for tests of strong vs. weak splitting that would be insensitive to corrections toq. Section V discusses why we anticipate that, for QCD, this particular test will be only partly successful at eliminating sensitivity toq eff . Finally, we offer our conclusions in section VI.

A. Simplifying assumptions
For simplicity, both for the sake of calculation and the sake of constructing as simple a thought experiment as possible to test weakly-coupled vs. strongly-coupled splitting, we will treat the underlying value ofq as constant in space and time. In contrast, the properties of quark-gluon plasmas produced in actual heavy-ion collisions depend on both.
In quoting results for large-N f QED, we will make a further simplification. This theory does not have the NLO double-log enhancements (1.10), but that does not mean that the value ofq relevant to splitting calculations is independent of formation length and so of energy. For QED, there is already (single) logarithmic dependence in the value of leadingorderq (0) . (For a very brief review qualitatively comparing and contrasting the QED and QCD energy dependence ofq in the language of this paper, see, for example, Appendix C of ref. [13].) Since our interest is not so much in the details of large-N f QED but in using it as an example for the calculation of overlap corrections, we will bypass this particular complication and simply treatq (0) as independent of energy. One could do better if interested: For QED (unlike QCD) the tentative assumption N f α EM (qE) 1/4 ) ≪ 1 of weak splitting in our calculation at large energy means that the QED plasma itself is also weakly coupled, N f α EM (gT ) ≪ 1. So, if desired, one could use a perturbative calculation of the relevantq (0) and its scale dependence. Since this is unnecessarily specific to the QED case, we've chosen to avoid complicating our presentation and will treatq (0) as constant.

B. More assumptions for the analysis in the main text
In the main text, we will focus on charge deposition rather than energy deposition because it is simpler to discuss at NLO. We leave discussion of energy deposition to appendix A.
We will also focus in the main text on the simplest derivations that will get us results for our example of large-N f QED, and leave discussion of more complicated situations to the appendices. In particular, we focus here on situations where the net "charge" carried by the parent of each splitting is unambiguously carried by just one, identifiable daughter of that splitting, with more general discussions in appendix A. Graphically, this corresponds  to simply following the original fermion line through the entire shower, such as the red line depicted in fig. 10, and then asking for the probability distribution for where that red line comes to a stop in the medium. If the splittings do not overlap (e.g. as in a calculation that used only leading-order results for splitting processes), there would be no problem: we could unambiguously identify the red line of fig. 10. At next-to-leading order, however, overlapping formation times allow the possibility of non-negligible interference between different finalstate fermions, as shown by fig. 11a. However, this interference is suppressed in large N. For both large-N f QED and large-N c QCD, it is suppressed (after summing over final fermion flavors or colors) compared to fig. 11b, for which the direct heir of the original fermion is unambiguous: the chance that the original fermion and the photon/gluon-produced pair have the same flavor or color is 1/N f or 1/N c suppressed, and distinguishable fermions cannot be confused with one another. In large-N c QCD, both figs. 11(a,b) are suppressed compared to fig. 11c after summing over final colors, and fig. 11c has no ambiguity regarding which daughter carries the quark number of the initial quark.
The advantage of cases where we can make this assumption is that, freed from having to follow anything but the single red line through any given shower like fig. 10, all that is needed to determine charge stopping can be packaged into a single differential rate dΓ/dξ for that red-line particle to reduce its longitudinal momentum by a factor of ξ. In a purely leading-order analysis, this rate would represent the bremsstrahlung process of fig. 12a. For a NLO analysis, it would also include the corrections to that rate due to overlapping formation time effects with the next splitting, such as depicted in figs. 12b and c. 6 The 12: Examples of processes which would contribute to the calculation of the differential rate dΓ/dξ at which a particle carrying a conserved charge (such as fermion number) degrades in energy: (a) leading-order bremsstrahlung; (b) the NLO effect of a subsequent splitting that has an overlapping formation time with the bremsstrahlung; (c) a virtual correction to the leadingorder bremsstrahlung (with overlapping formation time). In the case of large-N c QCD, the initial particle would be a quark (for study of charge deposition), but the particles in the real or virtual pair produced in (b) or (c) would be gluons.
rate dΓ/dξ has the added convenience that many of the contributions to overlapping double splitting cancel each other in dΓ/dξ, which simplifies computational work and is one reason we focus on charge deposition instead of energy deposition in this paper. We'll discuss the cancellations more concretely in section III. Here, we will just assume that the formula for dΓ/dξ is known through next-to-leading order.
The ambiguity of following the red line in fig. 11a is also the reason that the fermion energy loss rate dE/dz is ambiguous (beyond leading order) outside of the large-N limit.
In what follows, we first discuss a general equation (given the above assumptions) whose solution determines the statistical distribution of charge ρ(z) in terms of splitting rates dΓ/dξ. After, we will discuss the fairly simple solution for the moments of ρ if one ignores logarithmic factors in rates. We then present a trick for extending that analysis to the case where next-to-leading-order splitting rates contain terms involving single logarithms of energy, which will cover the case of large-N f QED. Finally, we mention (and then defer to the appendices) a different and computationally more difficult method that could be used in the case of more-complicated dependence on energy, such as double logarithms.
C. Integro-differential equation for the charge distribution ρ(z) Let ρ(E, z) be the distribution of charge deposition from an initial charge of energy E. To find an equation for ρ in terms of splitting rates, consider ρ(E, z+∆z) where ∆z is tiny, and think of traveling the distance z+∆z as first traveling ∆z followed by traveling distance z. In the first ∆z, the particle has a chance of not splitting at all, in which case the chance of traveling the remaining distance z will just be ρ(E, z). Alternatively, there is a chance that the particle splits in the first ∆z, reducing its energy to ξE, in which case the chance of traveling the remaining distance is ρ(ξE, z). In formulas, in section I.A of ref. [19].
where Γ(E) is the splitting rate for a particle with energy E. Re-arranging terms and taking the limit ∆z → 0, which can be rewritten as 3) The total LPM splitting rate Γ = (dΓ/dξ)dξ is typically infrared (IR) divergent from soft bremsstrahlung (at leading order in QCD and from NLO overlap corrections in large-N f QED), corresponding to ξ → 1 above. However, the factor ρ(E, z) − ρ(ξE, z) in the integrand of (2.3) will give an extra suppression as ξ → 1 that makes that ξ integral finite. So, even though we introduced IR-divergent quantities in the derivation of (2.3), the final equation is IR-safe.
The above argument for the ρ equation, as well as what we will do next, was inspired by related arguments made directly for the first moment ℓ stop of ρ in ref. [16]. It is also reminiscent of evolution equations explored in refs. [21,22] for the distribution N(E, x, t) of shower particles in longitudinal momentum fraction x as a function of time t=z, except that here we focus directly on the spatial distribution ρ of deposition and need not keep track of the detailed distribution N in x. 7 D. Warm-up: charge stopping with no logarithms

Results for moments
Our simple parametric estimate (1.5) showed that splitting rates scale with energy as E −1/2 . If we ignore logarithmic corrections (e.g. from energy dependence of α andq), the simple energy scaling allows us to to derive simple results for stopping lengths ℓ stop , and for other moments of the deposited charge distribution, in terms of dΓ/dξ. To see this, pull out the explicit energy dependence of the splitting rate by writing The only reason we introduce the symbol ν for the power 1 2 is that later it will help in generalizing this analysis to the case of single logarithms.) Because the rate scales as E −ν , the distance that the charge travels before stopping in the medium will scale as E +ν . Let's factor this scaling out of the probability distribution ρ(E, z) by replacing it with a probability distributionρ(z) in terms ofz ≡ z/E ν , normalized so that ρ dz =ρ dz. That is, (2.6) Plugging (2.4) and (2.6) into the ρ equation (2.3) gives which can be rewritten as As an aside, we note that numerical solution of this equation for leading-order QED splitting was used to generate the precise shape of ρ(z) shown as an example in fig. 7 (see appendix B). Now multiply both sides byz n and integrate overz to get (after integration by parts on the left-hand side of the equation) a recursive relationship between the moments of the probability distributionρ: giving z n = n z n−1 . (2.10a) From now on, we will assume that the charge q 0 of the particle initiating the shower is normalized to be q 0 =1, so that ρ is normalized such that 1 = 1. The case n=1 of (2.10a) then gives the formula in agreement with ref. [16]. 8 But the recursion relation (2.10a) derived here lets us find higher moments like z 2 and so σ = E ν z 2 − z 2 .

Measures of size of NLO correction
We'll now expand (2.10) to first order in α to get formulas for the measures ∆ℓ stop /ℓ (0) stop and χα of the relative size of corrections to showering from overlapping formation times.
To this end, decompose the splitting rate dΓ/dξ into its leading order piece (the usual LPM calculation for single splitting) and its NLO correction: The NLO piece is suppressed by O(α) compared to the leading-order rate, and the ≃ above just indicates that we are ignoring yet-higher order terms, suppressed by O(α 2 ) or more. Plugging (2.11) into (2.10b) and expanding ℓ stop ≃ ℓ (2.12) We find the structure of later equations to be a little clearer if we use the alternative notation where we define leading and next-to-leading order rate-weighted averages of any function f (ξ) as (2.14) Similarly, plugging (2.11) into the n=2 case A feature of (2.17) worth noting is that if the next-to-leading order rate dΓ (NLO) /dξ contains a piece that is any ξ-independent constant times the leading-order rate dΓ (0) /dξ, that piece will completely cancel between the two terms of (2.17) and so will not contribute at all to χα. So, for instance, if the NLO rate were to contain contributions corresponding to ξ-independent NLO corrections to theq appearing in the leading-order rate, such corrections toq would not contribute to χα. This was part of our earlier motivation to consider the ratio σ/ℓ stop and thence χα.

E. Charge stopping with single logarithms
We now generalize the previous analysis to a case that will cover QED in the approximations used in this paper. Consider the parametric formula (1.5) for the distance l rad between nearly-democratic splittings. Taking its inverse, the correspond splitting rate is In general, the leading-order splitting rate (dΓ/dξ) (0) is proportional to the α = α(µ) associated with the splitting vertex, like in the parametric formula above. As previously discussed, the renormalization and running of α means that logarithms ln(µ/Q ⊥ ) must therefore appear at next-to-leading order, similar to (1.16). Specifically, for large-N f QED, ref. [13] found by explicit calculation that the rate (2.11) has leading-order behavior dΓ (0) /dξ ∝ E −1/2 , like (2.4) and (2.18), but NLO corrections of the form 9 where here is the coefficient in the 1-loop renormalization group β-function for α EM . We'll find it convenient to rewrite (2.19) as In this section, we consider generally cases like (2.21) where the NLO corrections to the leading-order rate may have, in addition to terms with energy dependence E −1/2 , a term with energy dependence E −1/2 ln E proportional to the leading-order rate. Similar to (2.4), we will scale out all the explicit factors of E −1/2 by writing where NLO ′ means all the NLO terms except the ln E term. In principle, the coefficient " 1 4 β 0 " above could represent any coefficient of any sort of NLO single-log correction to the leading-order result, but we will use the notation 1 4 β 0 relevant to the QED application (2.21). ) − ln f (ξ) for any ξ-dependence f (ξ), the ξ dependence of the scale can be absorbed into the last term of (2.19). Here, the adjustments we need to make to our previous derivation of formulas for moments of the stopping distribution will only revolve around the dependence of the argument of the logarithm on E, not on ξ.
We can now use the following trick. Through next-to-leading order, This allows us to use the same formulas (2.10) that we derived for moments in the previous section, except that ν is now (2.24). Accounting for this difference in ν when we expand in α, (2.13) and (2.17) are modified to These are the formulas that we will use later with large-N f QED results for dΓ/dξ in order to produce figs. 6 and 8. For cases (such as QCD) where the NLO contribution to dΓ/dξ contains terms with more complicated energy dependence that E −1/2 ln E, we have not figured out a way to get formulas as simple as (2.25) and (2.26) for our tests of overlap corrections to shower development. One could turn to direct numerical simulation of the original equation (2.3) for ρ(E, z), but in the general case we have not so far figured out any particularly efficient way (numerically or otherwise) to solve the equation and cleanly isolate the expansion of that solution to first order in α. We present a different approach in Appendix C, which has been relegated to an appendix because (i) it is more complicated to derive and will be more complicated to implement numerically, and (ii) the previous results (2.25) and (2.26) are all that we will need for the large-N f QED results presented in this paper. Also, since we do not yet have full QCD results for dΓ/dξ, we leave the task of numerically implementing the procedure described in Appendix C to later work.
Note that, because all our measures of the size of NLO corrections are linear in dΓ (NLO) /dξ, one could use different techniques for different contributions to dΓ (NLO) /dξ and then sum together the results. So, if desired, (2.25) and (2.26) could be used to calculate the contribution to ∆ℓ stop /ℓ stop and χα from all terms of dΓ (NLO) /dξ except the ln 2 E term, and then just the ln 2 E contribution could be evaluated using some other method like Appendix C.
FIG. 13: Leading-order time-ordered diagram for e → γe. Blue represents a contribution to the amplitude and red represents a contribution to the conjugate amplitude. Repeated interactions with the medium are present but not explicitly shown. This diagram should be added to its complex conjugate by taking 2 Re[· · · ].
FIG. 14: Time-ordered interference diagrams for e → eēe in large-N f QED [13]. Complex conjugates should also be included by taking 2 Re[· · · ] of the above.

A. The relevant diagrams
In preceding sections, we have talked abstractly about the rate dΓ/dξ relevant to following the original charged particle through the medium in cases (like large N) where the original particle's heir is always identifiable, even at next-to-leading order. Here we will relate this dΓ/dξ, for large-N f QED, to the specific diagrams and results of ref. [13]. In the conventions of ref. [13], the leading-order rate dΓ (0) /dξ is given by the time-ordered interference diagram fig. 13 (plus its complex conjugate) for e → γe; the NLO corrections due to overlapping real-double splitting e → γe → eēe are given by fig. 14 (plus complex conjugates); and the NLO virtual corrections to single splitting are given by fig. 15 (plus complex conjugates). Each diagram represents a contribution to the rate for the process, which is the product of an amplitude (represented by the blue part) times a conjugate amplitude (represented by the red part). In each diagram, only the high-energy particles are explicitly shown, but each high-energy particle is interacting an arbitrary number of times with the medium, and a medium average is performed. See refs. [13,19,24] for more detail. The vertical photon lines with bars across them represent photons with longitudinal polarization in light-cone gauge, which produce instantaneous interactions in light-cone time, whereas the unbarred photon lines represent transverse photons.
All together, the dΓ/dξ discussed earlier in this paper is In the last equation, the two terms correspond respectively to figs. 14 and 15. ξ is called x e in ref. [13] and is the factor by which the original charge loses energy via single or overlapping double splitting. y e represents the momentum fraction of the pair-produced electron relative to the photon that produced it. The ∆ in ∆ dΓ indicates that one should subtract from fig. 14 for e → eēe the result one would get instead if one simply combined the leading-order formulas for the rate of e → γe with that for γ →ēe. In this way, double counting is avoided when using dΓ/dξ to generate shower development, and ∆ dΓ represents the corrections to the rate due to overlapping formation times. See section I.A of ref. [19] for more explanation. The ∆ dΓ in the third term of (3.1) indicates a similar subtraction [13] for fig. 15, where the second splitting is virtual.
Ref. [13] shows that there is a simple relationship (called there a back-end transformation) that relates the double-splitting diagrams (a,b,c,e,g) to the virtual diagrams (h,i,j,l,n) respectively. For the particular combination (3.1b) [in which the double-splitting diagrams are simply integrated over y e without any other factors], these diagrams in fact cancel each other and so do not need to be evaluated. Another relationship was found between (m) and (e), in a way that we have since discovered is related to (f) as 10 The "∆" designation in ∆ dΓ has been dropped because the diagrams where that distinction was important (diagrams a-c [19] and the related virtual diagrams) have canceled. In consequence, for calculations of the charge stopping distance in large-N f QED, we need evaluate just three of the NLO interference diagrams: the double-emission diagrams (d) and (f) and the virtual correction diagram (k).

B. Formulas and Numerical Integrals
To give a sense of the structure of the formulas for the relevant diagrams, and to help identify which formulas from ref. [13] are relevant, here are the top-level formulas from ref. [13], using that reference's notation. In particular, we write x e here instead of ξ, and is the momentum fraction of the pair-produced electron relative to the initial fermion of the double-splitting or virtual process, rather than relative to the photon. The leading-order rate is 11  11 The leading-order rate formula (3.5) is equivalent to the 1956 result by Migdal [15] if the latter is packaged into more modern, general notation. (See, for example, appendix C.4 of ref. [13] for a translation.) The QCD analogs, not shown here, trace back to BDMPS [2, 3] and Zakharov [4]. We refer the reader to ref. [13] for lower-level details. In particular, the symbols (ᾱ,β,γ, ζ) above are functions of (x e , y e ) related to splitting vertices and DGLAP splitting functions.
(The symbolᾱ should not be confused with the coupling α EM .) The various forms of the symbols (X, Y, Z, I, M, Ω) are functions of the relevant (x 1 , x 2 , x 3 , x 4 ) and arise in the evaluation of medium-averaged propagation of the high-energy particles between the splitting vertices in the diagrams. Renormalization of diagram (k) has been carried out in the MS-bar scheme (MS), and so α EM = α EM (µ) is the MS-bar renormalized coupling. The important thing to take away from the above formulas is that (3.6b) and (3.7b) involve performing a time-integral (∆t) of a complicated integrand. This integral must be done numerically because we have failed to do it analytically. Subsequently, an integral over y e or equivalently y e must be performed in (3.1b) and (3.7a) to obtain dΓ/dξ. Finally, integrals over ξ must be done to test NLO effects on shower development, such as (2.25) and (2.26). All told, that's a triple numerical integral (ξ, y e , ∆t) [equivalently (x e , y e , ∆t)] of a complicated integrand.
We found even the initial integration over ∆t to be numerically expensive for (3.7b). Presumably, more efficient numerical methods could be developed, but we took the following approach. We evaluated (3.7b) for an appropriately chosen mesh of points in the (ξ, y e ) plane. The function has integrable divergences in various limits such as ξ→0 or 1 and/or y e →0. Accounting for these divergences, we then figured out a way to numerically construct an interpolating function throughout the (ξ, y e ) integration region. Using that interpolation, we then performed the final integration over (ξ, y e ) by brute force with standard integration software (Mathematica [25]). Both the divergent limiting behavior and our method of interpolation are described in appendix D.

C. Results
Our main numerical results have already been displayed in figs. 6 and 8, which were generated by using the above formulas for dΓ/dξ in (2.25) and (2.26) respectively. The additional plot fig. 9 of the size of NLO corrections to dE/dz (which is well defined for large-N f QED) is given by A physical quantity, like ℓ stop , should not depend on the choice of renormalization scale µ if computed to all orders and expressed in terms of either (i) other physical quantities, or (ii) things that can in principle be extracted from other physical quantities. An example of the latter is MS-bar α EM (µ 0 ) at some specific choice of scale µ 0 that is expressed in terms of physical quantities, like our µ 0 ≡ (qE) 1/4 chosen in fig. 6. But generally, when one truncates a small-coupling expansion at some finite order, the result does have µ dependence-the size of the variation with µ is formally of order the size of yet-higher-order corrections which have not been calculated, which in our case would be NNLO. But the evaluation of ℓ stop through NLO order (for constantq (0) ) is an exception. Our expression for ℓ (0) stop + ∆ℓ stop (but neither term individually) is actually µ-independent provided one correspondingly uses the 1-loop approximation to the renormalization group equation for α(µ).
This feature is a consequence of the leading-order result ℓ to the 1-loop renormalization group equation In order for the stopping length to be µ-invariant through NLO, and given that ℓ where A and B are expressions that depend onq (0) and E but not on µ or α(µ). But then (4.1) means that the right-hand side of (4.3) does not depend at all on the value of µ if one uses 1-loop renormalization group equations for α(µ). The reason that our results for ∆ℓ stop nonetheless depend on µ is because the leading-order result ℓ (0) stop = A/α(µ) depends on µ implicitly through α(µ). Note that in fig. 6, we have chosen to keep the physics constant as we vary µ by taking the horizontal axis to be N f α(µ 0 ) rather than N f α(µ). Specifically, using (4.1) and the µ-independence of (4.3), we can rewrite (4.4) as . Note here that even for N f α(µ 0 ) = 1.0, the Landau pole is still two orders of magnitude away from the physics scale µ 0 , and so our calculations are sensible. Another way of saying this is that for some sorts of physics, like the running of the coupling, N f α(µ 0 ) = 1.0 can be considered a somewhat small coupling, whereas for the question of whether or not overlapping formation time effects are negligible in a medium, we have seen from figs. 6 and 8 that it is a moderately large coupling.

V. CAVEATS FOR χα TEST OF STRONG VS. WEAK-COUPLED SPLITTING
In the introduction, we motivated our discussion of the relative NLO correction χα in the expansion σ/ℓ stop ≃ (σ/ℓ stop ) (0) (1 + χα) as a (partly) successful method for testing the size of NLO corrections that could not be absorbed intoq. Here we explain the caveat "partly," which applies to the situation of double-log corrections toq in QCD.
The idea was that corrections toq would cancel between the numerator and denominator in σ/ℓ stop . This is true of energy-independent corrections toq, but the situation is more subtle for corrections that depend logarithmically on energy.

A. Warm up / review for a single logarithm
It will be instructive to imagine a theory where the NLO corrections toq depended on a single logarithm of energy that, so that δq(E) ≈ Aαq (0) ln(E/qτ 2 0 ) (for illustration only) (5.1a) since dΓ (0) /dξ scales like q/E. In this paper, depending on context, the symbol E sometimes refers (as above) to the energy of the parent of some splitting or overlapping doublesplitting somewhere in the shower, and it sometimes instead refers to the energy of the particle that initiated the entire shower. To avoid confusion, let's call the latter E 0 here, and write E = yE 0 for the parent of the NLO splitting in (5.1b), where y is the longitudinal momentum fraction of that splitting's parent relative to the total energy of the shower. Then (5.1) may be rewritten as The potentially large first term in (5.2a) is a constant addition toq and so the first term in (5.2b) should cancel in the calculation of the ratio σ/ℓ stop . We've seen that explicitly in the discussion after (2.17). The second term in (5.2b) does not generate large logarithms because, as discussed in the introduction with regard to (1.6), the late-time evolution of the shower (small y above) does not contribute significantly to stopping distances. In formulas, the effect of the ln y term in (5.2b) would be precisely the Avg[· · · ln ξ] (0) terms in (2.26) [with the constant " 1 4 β 0 " there identified as the 1 2 A of (5.1b)]. Those Avg[· · · ln ξ] (0) terms represent convergent integrals over ξ and do not generate any logarithmic enhancement of the result for χα.

B. Double logarithms
Now consider in contrast a double-log correction toq, as in QCD: 3) The first term above, which is the largest, is a constant and so will not contribute to σ/ℓ stop and so will not affect χα. The last term will not give any large-log contribution to σ/ℓ stop because stopping distances are not sensitive to small y. But the middle term can be expected to give a single-log enhanced contribution to σ/ℓ stop and so to χα. We characterized our proposal of studying χα as only partly successful because it manages to avoid double-log corrections (as well as avoiding all non-log corrections toq), but there will be an issue with sub-leading logarithms in QCD. For QED, there are no double logs, and so this difficulty does not arise there.

VI. CONCLUSION
In this paper, we have proposed various simple characteristics of in-medium showers that can be used to test the size of overlapping formation time corrections, in an attempt to distinguish weakly vs. strongly-coupled splitting pictures of shower development. Part of our motivation was to find a measure insensitive to corrections that can be absorbed intoq. We found one for QED, but our proposal seems like it would run into difficulty at the level of sub-leading logarithms in QCD. This underscores the importance of understanding how to account for and absorb subleading logarithms in ongoing research on the calculation and structure of overlapping formation time corrections in QCD.
We used large-N f QED as a concrete example in this paper, finding, for example, that the effects of overlapping formation times are modest (the weak-splitting picture appears good) for N f α EM (Q ⊥ ) ≃ 0.2, but is roughly a 100% effect for N f α EM (Q ⊥ ) ≃ 1. One is of course tempted to rashly speculate about how this might or might not translate to the the QCD analog C A α s =N c α s of N f α EM , but that would be premature. For one thing, the infrared behavior of the LPM effect in QED and QCD is very different. The completed set of rates needed for a (large N c ) QCD NLO calculation of shower characteristics will hopefully be available in the future, but the issue of subleading logarithms will then need to be understood.
In large-N f QED, the results of fig. 6 indicate that the net effect of overlapping formation times is to reduce the stopping length (corresponding to an increase in dE/dz). It would be nice to have a simple physical picture for understanding this sign of the result.
One may combine terms to rewrite (A1) as One could also use the ξ ↔ 1−ξ symmetry associated with having two identical daughters to simplify to ∂ε(E, z) ∂z but, for the sake of later generalization, we will find it useful to stay with (A2). Now consider the case where rates may be taken to scale with energy as E −ν , as in sections II D and II E of the main text. Write analogous to (2.6). This yields as the energy-deposition analog of (2.8). Taking moments of both sides leads to the recursion relation z n = n z n−1 Compare to (2.10). Note that in this appendix (except for section A 2) the notation z n refers to moments of the energy deposition distribution ε, whereas elsewhere in this paper it refers to moments of the charge distribution ρ.
One could now expand in α to obtain formulas for NLO corrections similar to (2.13), (2.17), (2.25), and (2.26), but we will not go into that level of detail here.

b. 1→2 splitting: multiple species
Now allow for multiple species, such as quarks and gluons. Showers initiated by different species i will generally produce different distributions ε i (E, z). The corresponding generalization of (A1) and (A2) is Here, for j=k terms of the double sum over species (e.g. g→gg in QCD), the overall factor of 1 2 in front of the integral is again a final-state identical-particle symmetry factor. For the j =k terms of the sum, the overall factor of 1 2 is canceled by the two permutations of jk in the double sum. For example, for i an electron in QED, the sum over j and k on the right-hand side would contain two terms which would be exactly equal: e→γe and e→eγ. The summation in (A7) is just a compact way of simultaneously accounting for the cases of identical and non-identical daughters.
Scaling out the energy dependence and taking moments of the equation, one finds As indicated by the last line, it is convenient to write the relation in terms of a matrix M (n) acting on species space, with matrix elements The recursion relation for the moments is then given in terms of the matrix inverse: The n=1 case reproduces (for ν = 1 2 ) the energy stopping length formula derived in ref. [16] for the case of quarks and gluons.

c. 1→2 and 1→3 splitting
The point of this paper is to be able to account for next-to-leading order corrections to shower development, which are generated when two consecutive splittings overlap. As discussed in section I.A of ref. [19], overlap corrections to double splitting should be treated as effectively a type of 1→3 splitting in order to analyze shower development in the framework of classical probability theory. 12 In addition, there are also direct 1→3 splitting processes at NLO due, for example, to the 4-gluon vertex in QCD [26] or interactions involving intermediate longitudinally-polarized gauge bosons [13].
Let dΓ i→jk and dΓ i→jkl be the rates for 1→2 and effective 1→3 processes respectively. 12 Depending on whether overlap effects reduce or enhance the double-splitting process, the "rate" assigned to this effective 1→3 splitting process could be negative. But that does not cause any difficulty for the type of analysis considered in this paper.
The generalization of (A7) to include the 1→3 processes is where we have now used δ-functions to write the integrals over daughter momentum fractions ξ in a more symmetric form. For rates that scale with energy as E −ν , the recursion relation for the moments is again (A10) but now with

. Cancellation of QCD power-law divergences
In the above formulas, we have (i) grouped overlap corrections to actual double-splitting (e.g. fig. 14) into dΓ i→jkl and (ii) grouped the corresponding NLO virtual corrections to single splitting (e.g. fig. 15) together with leading-order single splitting into dΓ i→jk . However, in QCD, these LPM rates individually have severe power-law (not just logarithmic) infrared divergences that make the two terms in (A12) separately infinite. Specifically, refs. [19,24] show for large-N c QCD that 13 (See section I.D of ref. [19] for a qualitative discussion.) As noted there, this blows up fast enough as ξ 3 → 0 to have a divergent effect on energy loss if NLO virtual corrections to single splitting are ignored. In particular, (A15) means that the Avg 1→3 [· · · ] in (A12) has a power-law divergence from, for example, the ξ 3 → 0 region of integration in (A14). It is firmly expected that this QCD power-law divergence will cancel in the sum in (A12) of Avg 1→3 [· · · ] and the NLO piece of Avg 1→2 [· · · ], leaving only the known QCD double-log divergence. This is something to verify in the future when results are available for the NLO virtual corrections to QCD single splitting. (The double log term itself needs to be treated using a more sophisticated method for computing the moments z n , as discussed in section II F.) 2. Charge deposition without assuming large N Using the same methods as above, we can also generalize the charge deposition discussion of section II D to handle the case where N f is not large. Beyond leading-order, we must then follow all of the charged daughters of each splitting because of interference effects like fig.  11a.
To adapt (A11) to charge deposition, it's useful to normalize each ρ i (E, z) so that its integral (i.e. 1 i ) gives the relevant charge of species i. Then the charge deposition analog of (A11) is However, in QED (and analogously in QCD), charge conjugation invariance gives and so one may rewrite (A16) as an equation just for ρ e (E, z): For rates that scale as E −ν , the corresponding recursion relation for the moments of ρ e would be z n = n z n−1 In this appendix, we discuss how to numerically solve for the full leading-order chargestopping distribution ρ (0) (z), which is interesting but not necessary for anything else in this paper. We will not have to restrict attention to the large-N f limit of QED in this appendix: large N f was used to simplify the complexity of NLO calculations but this appendix only deals with leading order.

Numerical Procedure
In particular, we will solve the scaled equation (2.8), which is directly relevant only to leading-order calculations (given our assumptions) but not to NLO (which contains log dependence on energy). Note that (2.8) is linear inρ, and so solving that equation does not by itself determine the overall normalization ofρ. But we may numerically find a solution with any normalization and then, as a final step, compute the normalization N = ∞ 0 dzρ(z) and then rescaleρ(z) by a factor of 1/N to makeρ into a normalized probability distribution.
To solve numerically, discretize z and rewrite (2.8) as for small steps −∆z and ν = 1 2 . This version of the equation lets us take a solution for largẽ z (sayz > some Z) and extrapolate it step by step to smallerz. All we need to start the process is an approximate (un-normalized) solution to (2.8) in the large-z limit.

Large-z solution
We expect the charge deposition distribution should decay rapidly (e.g. exponentially) for large z, so let us make a WKB-inspired rewriting where we will treat the exponent W (z) as large. Plugging (B2) into (2.8), We've explicitly written the superscripts "(0)" at this point to remind us that we are only solving the equation for the leading-order result. We've made that reminder because it will be important here that the LPM suppression of leading-order soft bremsstrahlung rates is very different in QED and QCD, such that the total LPM bremsstrahlung rate Γ (0) is infrared divergent in QCD but not QED. For finite total Γ (0) , we assert (and will justify a posteriori) that the second term in the integrand of (B3) can be ignored whenz is large, leaving Putting this solution back into (B3), one finds that, forz ≫ 1/Γ (0) , the second term in the integrand is indeed ignorable compared to the first unless 1−ξ (Γ (0)z ) −1 . That leaves only a parametrically small portion of the ξ integration where the second term is important. As long as the integrals of the first and second terms are separately convergent, that means that the effect of the second term on the integral is negligible in the largez limit. For leading-order QCD, one needs a different analysis, which we will not present here. We have not bothered to go to higher order in our WKB-like expansion to determine whether there is anyz dependence to the pre-factor of the exponential behavior (B6). For large enoughz, an algebraic pre-factor will not vary significantly over the range 1/Γ (0) it takes for the exponential to drop drastically. So our numerics will not be very sensitive to that pre-factor provided we only use (B6) for large enoughz. This insensitivity may be checked by varying the cut-off Z chosen to switch between the asymptotic form forz > Z and the numerical evolution using (B1) forz < Z.

Result
Using the above procedure with the QED leading-order LPM bremsstrahlung rate (3.5) gave the particular curve shown in fig. 7. 14 Since this is a precise numeric result, we now show it in a little more detail in fig. 17.
Appendix C: One method for charge stopping with double logarithms, etc.

Basic equations
Instead of directly writing an equation (2.3) for the charge deposition, consider following the original charge through time, even before it is deposited, from energy E initially to some smaller xE at time t = z. Let P (x, E, t) be the probability distribution for x at time t. We can find an equation for P (x, E, t) somewhat similar to how we found the equation (2.3) for ρ(x, E), and very similar to the method used by ref. [21] for gluon distributions. The basic equation for the evolution of P over an infinitesimal time interval ∆t is (C1) The first term is the probability that nothing happened in the small time interval [t, t+∆t] multiplied by the then-unchanged P (x, E, t). The second term is the probability density P (y, E, t) that the particle had some longitudinal momentum fraction y > x at time t, convolved with the probability that a splitting took it from y to x during the small time interval [t, t+∆t]. Doing the y integral, and re-arranging terms and taking the limit ∆t → 0, This can be rewritten in the infrared-safe form where, in order to compactly combine integrals, we have adopted the physically-appropriate convention that P (x, E, t) ≡ 0 for x > 1.
The initial condition on the time-evolution (C3) is that the particle start with energy E: As noted in the similar context of ref. [21], the distribution P (x, E, t) will have a δ-function piece corresponding to stopped particles: where P stopped (E, t) in our case is the chance that the original charged particle has already stopped by time t. One can use this to recover the charge deposition distribution ρ(x, z) that will be needed to calculate, for example, σ/ℓ stop . The total probability that the charge is still moving at time t is where 0 + represents a positive infinitesimal. The charge deposition distribution ρ(E, z) [where t and z are interchangeable in our discussion] is then The moments z n of that distribution are then (integrating once by parts) (C9)

Leading-Order Equation
As in the main text, assume that the leading-order splitting rate scales with energy as E −1/2 without any logs. Even though we then already know how to directly calculate the leading-order values of the moments z n using the method of section II D, we will still, it turns out, need to know the leading-order solution P (0) (x, E, t) to (C3) in order to find the NLO correction P (NLO) (x, E, t) to P , which we will need to find the NLO correction to moments. So let's first discuss how to find that leading-order solution. Similar to refs. [21,22], we can use the E −1/2 scaling of the leading-order splitting rate (and therefore E 1/2 scaling of distance and time scales) to simplify our basic equation (C3) by writing (C10) with initial conditionP (0) (x, 0) = δ(1 − x). We have not found any closed form solutions, but (C11) could be solved numerically. We will assume one has such a solution in hand and now discuss how to obtain from it the NLO correction.
3. Finding P (NLO) from P (0) 15 Our (C11) is similar in form to eq. (4) of ref. [21]. One difference is that their analysis accounts for all daughters of each splitting g → gg of a gluon shower because they are interested in the number distribution of gluons. We would also need to do this to compute what we call "energy" deposition by this method. But here we only need to track the fate of the single red-line particle in fig. 10 to compute charge deposition ρ(z, E) (at least in large-N theories, as discussed in the main text). Additionally, there are some inessential normalization differences associated with how we define dΓ from dΓ [related to their K(z) from dP br /dz dτ ] and how we define ourt from t [related to their τ from t].
Making no assumptions about the energy dependence of the NLO piece dΓ (NLO) /dξ of the splitting rate, we now show how to determine P (NLO) from P (0) . Rather than returning to the general equation (C3) for P (x, E, t), we find it easier to find P NLO by a direct probability argument. First, as a matter of language, let's call dΓ (0) /dξ the rate for leading-order splittings and dΓ (NLO) /dξ the rate for "NLO splittings" (e.g. the correction from overlapping double splittings, or the virtual correction to single splittings). Then imagine a shower composed of any number of leading-order splittings plus zero or one or more NLO splittings, and correspondingly decompose This is not quite the same as (C12) because P (0) is the probability if the only type of splitting possible were leading-order splittings. P (0 NLO splits) , in contrast, imagines a world in which NLO splittings are also possible but, by random chance, none have occurred in the time interval [0, t]. It differs from P (0) by having to account for the probability that no NLO splittings occurred. Before examining the difference between P (0) and P (0 NLO splits) in more detail, it will be helpful to first figure out P (1 NLO split) , depicted in fig. 18. In formulas, this picture translates to The first factor in the integrand gives the probability of making it to time t ′ with only leadingorder splittings which take the energy of the original particle down from E to yE. The second factor gives the probability that there is then a NLO splitting in the interval [t ′ , t ′ +dt ′ ] that takes the energy from yE to ξyE. The final factor (in square brackets) is related to the probability of then making it from there to the later time t with only leading-order splittings taking the energy from ξyE to the final value xE = (x/ξy) × (ξyE). The reason for the normalization factor (ξy) −1 inside the square brackets is because the left-hand side of (C14) represent a probability distribution for x whereas the P (0 NLO splits) ( x ξy , ξyE, t − t ′ ) on the right-hand side represents a probability distribution for x/ξy, and so one must convert.
In (C14), the explicit factor dΓ (NLO) /dξ is already next-to-leading order in α. Since we are ultimately only interested in the expansion (C12) of probabilities characterizing the shower to next-to-leading-order, we may therefore approximate the other factors P (0 NLO splits) in (C14) by P (0) . So Now let's return to the first term of (C13), P (0 NLO splits) , and its difference from P (0) at first order in dΓ (NLO) . To evaluate this difference, consider fig. 19, which is the analog of fig. 18 but with no NLO splitting occurring during the time interval [t ′ , t ′ +dt ′ ]. The where, in order to combine integrals, we have again adopted the convention (C4). 4. NLO corrections to moments z n Similar to our notation ℓ stop ≃ ℓ (0) stop + ∆ℓ stop , we'll write the expansion of other moments of the charge deposition distribution ρ(z, E) to NLO as z n ≃ z n (0) + ∆ z n . (C19) Putting (C18) into (C9), Because dΓ (0) and P (0) have simple scaling with E, we will find it convenient to scale out the same factors of E 1/2 from the NLO rate as in the main text: Using (C10), eq. (C20) above can be written ∆ z n = −n The constraint (C4) allows us to replace the lower limit of integration on the y integral by zero. Continuing to use (C4), the x integration in (C25) can then be accomplished as follows: at leading-log order, where (also at leading-log order) one may use parametric estimates for the formation times in the argument of the logarithm. We find numerically that this same formula works for the NLO single-splitting correction (D1). Moreover, we find it works whenever the formation time t form,ye for pair production from the photon is parametrically small compared to the formation time t form,xe for the bremsstrahlung, which includes the limit y e → 0 as well as x e → 1. Using numerics, we have also extracted the constant under the logarithm for this limit, finding that f (x e , y e ) approaches For the x e →0 divergence, we do not know of a physical argument like the one of ref. [13] for the previous limit (D3), and we have not tried to derive this limit analytically from the general formula for f (x e , y e ). However, we have found by numerical experimentation that f (x e , y e ) approaches F 0 (x e , y e ) ≡ −N f α 2 EM c 0 P γ→e (y e ) for x e →0 with y e fixed, where c 0 is a constant whose numerical value is approximately .
Next, we combined the limits F 0 and F 1 to make a positive weighting function F that correctly approximates the magnitude of f in all of the divergent limiting cases. The goal was to create a new function g(x e , y e ) ≡ f (x e , y e ) F (x e , y e ) (D10) that would not be divergent. By trial and error, trying to choose a g that looked reasonable, we settled on the choice F (x e , y e ) ≡ [F 1 (x e , y e )] 2 − F 1 (x e , y e ) F 0 (x e , y e ) + [F 0 (x e , y e )] 2 .
The resulting finite function g(x e , y e ) is shown in fig. 20b. Note that it has an unfortunate finite directional singularity for (x e , y e ) → (0, 0), where the value of g depends on how one approaches the limit. One can trace this back to a similar directional singularity in how f diverged as (x e , y e ) → (0, 0) in fig. 20a. The directional singularity in fig. 20b complicates good numerical interpolation of g from a finite mesh of values. We therefore looked for a y e -dependent change of the x e variable, x e →x e (y e ), that would map the interval [0, 1] into itself while shifting the location of the ridge in fig. 20b to avoid running into the corner. To map [0, 1] into itself, we considered a transformation of the formx e (y e ) = 1 − x 0 (y e ) x e x 0 (y e ) + 1 − 2x 0 (y e ) x e , which maps x e = x 0 (y e ) tox e = 1 2 . By experimentation, we chose  Fig. 20c then shows g as a function of (x e , y e ) instead of (x e , y e ). This function is smooth enough for numerical interpolation. Our procedure was to evaluate f for a mesh of points in (x e , y e ) space, from which we obtain g on that mesh of points, and then interpolate to get g in the entire (x e , y e ) region. The original function f (x e , y e ) can then be evaluated from this interpolation by inverting all of the steps described above.
One difference here is that the only divergence of f is for x e →1. However, there is √ y e or √ x e behavior in other limits, which is best to also remove before interpolation. So we again look at all the limiting cases. We find  x eq E 1/2 with to get fig. 21b. Finally, here we take (D12a) with x 0 (η) = 1.2 y e (1−y e ) (D18) to get fig. 21c. Any reader going through this appendix may be forgiven if they think it's all a bit convoluted. We agree and would be happy to have a more efficient procedure (in terms of our time or the computer's).