A local Monte Carlo framework for coherent QCD parton energy loss

Monte Carlo (MC) simulations are the standard tool for describing jet-like multi-particle final states. To apply them to the simulation of medium-modified jets in heavy ion collisions, a probabilistic implementation of medium-induced quantum interference effects is needed. Here, we analyze in detail how the quantum interference effects included in the BDMPS-Z formalism of medium-induced gluon radiation can be implemented in a quantitatively controlled, local probabilistic parton cascade. The resulting MC algorithm is formulated in terms of elastic and inelastic mean free paths, and it is by construction insensitive to the IR and UV divergences of the total elastic and inelastic cross sections that serve as its basic building blocks in the incoherent limit. Interference effects are implemented by reweighting gluon production histories as a function of the number of scattering centers that act within the gluon formation time. Unlike existing implementations based on gluon formation time, we find generic arguments for why a quantitative implementation of quantum interference cannot amount to a mere dead-time requirement for subsequent gluon production. We validate the proposed MC algorithm by comparing MC simulations with parametric dependencies and analytical results of the BDMPS-Z formalism. In particular, we show that the MC algorithm interpolates correctly between analytically known limiting cases for totally coherent and incoherent gluon production, and that it accounts quantitatively for the medium-induced gluon energy distribution and the resulting average parton energy loss. We also verify that the MC algorithm implements the transverse momentum broadening of the BDMPS-Z formalism. We finally discuss why the proposed MC algorithm provides a suitable starting point for going beyond the approximations of the BDMPS-Z formalism.


Introduction
Most generally, the notion 'jet quenching' is currently used to characterize a broad range of experimental observations including the modification of high-p T single inclusive hadron spectra, jet-like particle correlations and reconstructed jets in nucleus-nucleus collisions. Jet quenching was discovered at RHIC via measurements of single inclusive hadron spectra, and the phenomenon was characterized extensively on the level of two-particle near-side and back-to-back high-p T correlation functions and particle yields associated with trigger particles [1,2]. Two-particle correlations displaying very similar features were also seen at the same time in nucleus-nucleus collisions at the ten times lower center of mass energy of the CERN SPS [3], whereas the single inclusive hadron spectra at the CERN SPS do not show the dramatic suppression up to a factor 5 observed at collider energies [4]. In recent years, a strong effort has gone into studying jet quenching at the highest experimentally accessible transverse momenta where one may hope to establish the most direct link between the rich jet quenching phenomenology and a partonic explanation rooted in QCD.
In this context, we mention that first preliminary results on reconstructed jet measurements have become available at RHIC [5][6][7] within the last two years. With the much wider kinematic reach accessible at the LHC, numerous novel opportunities for studying jet quenching emerge now. Data from the first exploratory heavy ion run have shown already that the nuclear modification of charged hadron spectra is somewhat stronger than at RHIC and that it persists beyond p T = 20 GeV [8]. Soon, the kinematic range of these measurements will be extended by a large factor, and much more detailed information about quenching of high-p T particles and particle correlations will become available. Moreover, first measurements of reconstructed jets in heavy ion collisions at LHC indicate already that also jets of order 100 GeV display significant medium-modifications. In particular, samples of reconstructed dijets display an energy imbalance distribution that is much wider than in the absence of a nuclear environment [9,10]. These measurements indicate that the quenching of reconstructed jets is accompanied by the medium-induced production of many soft particles [11]. At present, our still incomplete theoretical understanding of jet quenching is largely based on the picture that highly energetic partons produced in dense QCD matter are degraded in energy due to elastic and inelastic interactions with the surrounding medium prior to hadronization outside the medium [12][13][14][15][16]. This picture is supported in particular by data on single inclusive hadron spectra and particle correlations. The coming years are likely to show a strong interplay of experimental and theoretical efforts to characterize jet quenching on the level of multi-particle final states and reconstructed jets with the aim of further constraining the microscopic dynamics of this phenomenon and drawing conclusions about the properties of the QCD matter by which it is induced.
Monte Carlo tools have well-recognized advantages for the phenomenological analysis of highp T multi-particle final states. They are the method of choice for formulating the evolution of a parton shower with minimal kinematic approximations and exact implementation of conservation laws. They are also best suited for interfacing this partonic evolution with the hadronic final state. Moreover, the fact that they generate not only event averages but also event distributions of final state particles meets an obvious experimental demand and allows for the interfacing with modern jet finding techniques [17]. To satisfy these experimental and theoretical needs for the study of heavy ion collisions, several Monte Carlo tools for the simulation of jet quenching have been developed in recent years. Some of the available tools are full event generators that supplement standard 'vacuum' final state parton showers with models of medium-induced gluon radiation tailored to analytical calculations of medium-induced parton energy loss. Hijing [18,19], Q-Pythia [20], Q-Herwig [21] and Pyquen/Hydjet++ [22,23] fall into this class. Other approaches modify the Pythia parton shower, e.g., to implement the picture of a medium-modified Q 2 -evolution as in YaJEM [24,25], or to implement rate equations based on a perturbative calculation of partonic energy loss as in Martini [26]. Finally, Jewel [27] aims at formulating a stand-alone final state parton shower that interpolates between three analytically known limits, namely the vacuum parton shower in the absence of medium effects, the analytically known limit of energy loss via elastic multiple scattering, and radiative energy loss. In its current version, however, radiative energy loss is modeled similar to other efforts ad hoc in terms of medium-modified splitting functions. A more detailed discussion of the current status of MC tools for jet quenching can be found in Ref. [28].
The 'vacuum' parton showers used in MC event generators like Pythia [29], Herwig [30] and Sherpa [31] are faithful representations of the theory of Quantum Chromodynamics (QCD). They resum to leading logarithmic accuracy the large logarithms associated with collinear gluon emission, and they thus implement with known accuracy and without additional model-dependent input analytically known features of QCD. In contrast, the MC tools for jet quenching listed above are phenomenological models. They may tailor some numerical steps according to QCD-based analytical calculations, but these QCD-based results do not define the MC tool up to controlled accuracy, they solely motivate physical choices in a more complex (and more complete) dynamical procedure. This is a perfectly legitimate approach that meets the demand of a broad range of applications. We argue, however, that it is also of interest to complement these pragmatic approaches with a conceptual exploration of whether a MC algorithm of jet quenching can be formulated as a faithful implementation of QCD-based calculations of parton energy loss. Establishing such a clearer connection between MC tools and analytical QCD-based knowledge of jet quenching may be important for constraining the fundamental QCD properties of matter that induce the observed jet quenching phenomena. Moreover, as we shall discuss in detail in section 7, such a faithful MC implementation provides a suitable starting point for overcoming many of the technical limitations of the state of the art of analytical parton energy loss calculations. With this motivation, we present in the present paper a MC tool that provides with controlled accuracy a local and probabilistic implementation of the BDMPS-Z formalism of medium-induced radiative parton energy loss.
The BDMPS-Z formalism [32][33][34][35] is historically one of the first QCD-based calculations of medium-induced radiative parton energy loss in the high energy limit. Its path-integral formulation that we recall in section 2, provides the generating function for formulations of radiative parton energy loss in terms of an opacity expansion [36,37]. Also, other formulations of medium-induced radiative parton energy loss [38,39] are known to display the same medium-dependencies as the BDMPS-Z formalism ( for a more complete overview, see Ref. [12]). In short, the most widely used radiative parton energy loss calculations are closely related to the BDMPS-Z formalism. Moreover, all existing analytical results, as well as generic physics reasoning, point to the dominant role of the so-called non-abelian Landau-Pomeranchuk-Migdal (LPM) effect in medium-induced gluon radiation, and this destructive quantum interference effect is accounted for in the BDMPS-Z formalism. We therefore expect that a MC implementation of the BDMPS-Z formalism can provide more general guidance as to how medium-effects should be formulated in a MC parton shower.
We note as an aside, that the BDMPS-Z formalism does not provide all the information that enters a final-state parton shower. For instance, the BDMPS-Z formalism has been derived for a relatively limited kinematic range only (see discussion in section 2), and it does not specify whether and how the angular ordering prescription of a vacuum parton shower should be changed in the medium. For recent work on this latter question, see Ref. [40,41]. The present paper will not address these advanced issues. To the extent to which future studies of radiative parton energy loss result in improvements of the BDMPS-Z formalism, it will be interesting to explore whether these refinements can be incorporated in modifications of the MC algorithm discussed here.
A priori, it is unclear whether destructive quantum interference such as the non-abelian LPM effect can be recast in a local probabilistic MC implementation of controlled accuracy. A prominent example in which destructive quantum interference can be formulated indeed in terms of a probabilistic prescription is the angular ordering condition of the vacuum parton shower. In general, however, quantum interference effects need not be in one-to-one correspondence with a local and probabilistic procedure. In a previous paper, we had pointed out [42] that the concept of formation time can be identified unambiguously in the BDMPS-Z formalism and that it could play the same role for the probabilistic implementation of medium-induced quantum interference as does angular ordering for implementing destructive interference of gluon production processes in the vacuum. In the present paper, this basic idea is worked out in full technical detail. It will also become clear why some elements of our original proposal have to be modified to arrive at a faithful implementation of the BDMPS-Z formalism.
Our paper is organized as follows: We first identify the main building blocks of the proposed MC implementation by analyzing in section 2 the BDMPS-Z formalism in the opacity expansion. Based on this analysis, we discuss in section 3 a simplified MC algorithm that does not trace yet the kinematic dependences of parton splitting, but that accounts for the formal BDMPS-Z limits of totally coherent and incoherent gluon production on the level of total radiated particle yields. Section 4 discusses how this elementary algorithm extends naturally to a full MC implementation of the BDMPS-Z formalism. In sections 5 and 6, we demonstrate that the proposed MC algorithm provides indeed a quantitatively controlled implementation of the BDMPS-Z formalism. Finally, we discuss in the outlook of section 7 the perspectives for further uses and developments of this MC tool.

Time-scale for medium-induced interference in the opacity expansion
Medium-induced gluon radiation is expected to be the dominant energy loss mechanism of highly energetic partons in QCD matter. Several groups have calculated the corresponding mediuminduced gluon energy distribution ω dI dω in the kinematical regime [32][33][34][35][36][37][38] where the energy E of the projectile parent parton is much larger than the energy ω of the radiated gluon, which is much larger than its transverse momentum k and the transverse momentum transfers q i from scattering centers in the medium.
In this section, we recall first that to each order in opacity [36,37], the double differential medium-induced gluon distribution ω dI dω dk can be written in terms of two classes of elementary cross sections (called R and H and defined below), multiplied by weighting factors that interpolate between limits of coherent and incoherent particle production. We emphasize that the scales of interpolation between coherent and incoherent particle production are set by inverse transverse energies that have an interpretation as formation times. They will play a central role in the algorithm proposed in section 4.

Medium-induced gluon radiation in the high energy limit
Our aim is to specify a Monte Carlo algorithm that implements the double differential mediuminduced gluon energy distribution ω dI dω dk , derived first by Baier, Dokshitzer, Mueller, Peigné and Schiff (BDMPS) [32,33] and independently by Zakharov [34,35] in the eikonal approximation (2.1). As a preparatory step, we summarize here information about ω dI dω dk that will be needed in the following discussion. For a medium of finite size, the distribution ω dI dω dk of radiated gluons can be written in the compact path integral formulation [36] Here, the right hand side of (2.2) contains several internal variables (u, y, r, y l ,ȳ l ), which do not relate directly to measurable quantities. The longitudinal coordinates y l ,ȳ l result from integrating over the ordered longitudinal gluon emission points in the amplitude and complex conjugate amplitude of a multiple scattering cross section. The two-dimensional transverse coordinates u, y and r emerge in the derivation of (2.2) as distances between the positions of projectile components in the amplitude and complex conjugate amplitude [36]. In the following, we discuss in more detail how the hard 'projectile' parton, the 'target' medium, and the interaction between both is accounted for by equation (2.2). Characterization of the medium: A partonic projectile that interacts perturbatively with the medium exchanges gluons with some components of the target. The momentum transfer between projectile and target can involve both transverse momentum q and longitudinal momentum q l . In radiative parton energy loss calculations based on the high-energy approximation (2.1), the transverse momentum transfer dominates, |q| q l . This motivates a description of the target in terms of a collection of colored static scattering potentials A(q), [43,44] eikonal q l q, This approximation neglects recoil effects, and thus it automatically neglects collisional energy loss. To treat collisional and radiative energy loss on the same level, one would have to undo the approximation (2.3).
In equation (2.2), the scattering potentials A(q) enter the gluon energy distribution in the form of the so-called dipole cross section Here, |A(q)| 2 characterizes the differential elastic cross section with which the projectile parton transfers a transverse momentum q to a single scattering center in the medium. In the gluon energy distribution (2.2), this quantity is always multiplied by the density n(ξ) of scattering centers along the trajectory of the projectile. For notational simplicity, we focus in the following on a homogeneous density distribution of scattering centers within a box of length L, that is n(ξ) = n 0 , for 0 < ξ < L , 0 , for ξ < 0 or L < ξ . (2.5) Our discussion generalizes to arbitrary density profiles, but we shall not provide details about this generalization in the present work.
Initializing the parent parton: The lower bound ξ 0 of the y l -integral of (2.2) denotes the time at which the high energy parton is produced. The parton is produced either at some finite time, which we set to ξ 0 = 0, or it is produced in the infinite past. These two initializations correspond to different physics scenarios: If the parton is produced in a hard interaction, then it is produced at a finite production time, which we set to ξ 0 = 0. Even in the absence of a medium, partons produced in a hard collision branch as a consequence of their virtuality. Equation (2.2) contains information about this vacuum splitting, since it leads in the absence of a medium to where C R = C F for a projectile quark and C R = C A for a projectile gluon. In (2.6), the notation N = 0 stands for the zeroth order in the opacity expansion, which corresponds to the case n(ξ) = 0, where medium effects vanish. The result (2.6) can be identified with the LO g → g g and q → q g vacuum splitting functions in the form that these splitting functions take in the eikonal limit (2.1).
• ξ 0 = −∞ The condition ξ 0 = −∞ initializes a parton that has propagated for an infinitely long time without branching, prior to possibly interacting with the medium for times ξ ≥ 0. In the absence of a medium, this parton will never branch, In this sense, the parent parton propagates as if it were 'on-shell'. Because of confinement, a colored parton does not propagate forever and this situation will never be realized in a physical process in the vacuum. But it is a relevant limiting case for understanding the physics contained in (2.2) Characteristic interaction terms: In the following subsections, we shall demonstrate that the terms related to vacuum radiation and medium-induced radiation can be identified unambiguously in the radiated gluon energy distribution (2.2) even outside the incoherent limit. In preparation for this analysis, we here define the kinematic dependencies which signal vacuum radiation and medium-induced radiation.
Perturbative splittings in the vacuum result in a characteristic 1 k 2 -distribution of the daughter gluons, with the transverse momentum measured with respect to the direction of the high energy parent parton. As vacuum radiation term, we shall identify the term which appears for instance in equation (2.6). Consistent with vacuum radiation, this term does not depend on medium properties. If a gluon, produced by vacuum radiation, scatters incoherently on N scattering centers which transfer transverse momenta q i respectively, then the transverse momentum distribution of the gluon will be shifted to We will refer also to terms of the form (2.9) as (shifted) vacuum radiation.
In the eikonal approximation (2.1), the basic cross section for medium-induced gluon radiation in potential scattering with momentum transfer q between target and projectile can be written as Here, |A(q)| 2 characterizes the differential elastic scattering cross section with which the projectile parton interacts with the static potential, and R(k; q) is the Bertsch-Gunion term [45] R(k; q) = q 2 k 2 (k + q) 2 , (2.11) which denotes the distribution of gluons of transverse momentum k, produced in a single incoherent interaction of a high energy parton with a colored scattering potential transferring transverse momentum q. The Bertsch-Gunion term characterizes medium-induced radiation. Consistent with this notation, R vanishes in the absence of medium effects, that is for q = 0. If a gluon, after being produced incoherently on one scattering center, scatters subsequently incoherently on N − 1 other scattering centers, then the Bertsch-Gunion momentum distribution is shifted to This is the incoherent (i.e. probabilistic) result of multiple elastic scattering.
In analyzing the gluon energy distribution (2.2), we shall also encounter medium-induced radiation terms of the form (2.13) These terms result when N scattering centers act coherently in a single gluon production process. They will be found in regions of phase space, where the formation time of the gluon is too long to resolve the N scattering centers individually. Of course, radiation terms in which gluons are produced in the coherent scattering on N scattering centers prior to rescattering incoherently on M other scattering centers can be found also. These terms are of the form (2.14) In the following subsections, we analyze the gluon energy distribution (2.2) in the opacity expansion. In doing so, we substantiate the table In particular, we specify how interference effects interpolate between incoherent elementary processes of the form H and R. This will be the basis for proposing a MC algorithm. terms of form R and H Table 1: Summary of characteristics of the ξ 0 = 0 and ξ 0 = −∞ cases 2.2 Interference effects for medium-induced gluon radiation (case ξ 0 = −∞) As discussed above, the gluon energy distribution (2.2) for a parton initialized at time ξ 0 = −∞ allows us to study the interference of different sources of medium-induced radiation in a limiting case, in which complications due to vacuum radiation are absent.
The following analysis relies on the opacity expansion. This is an expansion of the integrand of (2.2) in powers of the density of scattering centers n(ξ) times the effective scattering strength σ(r) of a single scattering center. The opacity expansion amounts to an expansion in powers of dξ n(ξ) V tot = n 0 L V tot , where V tot characterizes the cross section presented by the scattering potential A(q), (2. 15) In practice, the opacity expansion of (2.2) results in integrations over the transverse momenta q 1 , ...,q N , which are weighted by the differential elastic scattering cross sections |A(q 1 )| 2 , ..., |A(q N )| 2 , but which do automatically factorize into powers of V tot . For this reason, the N -th order of opacity is obtained most easily by collecting all terms of order (n 0 L) N .

N = 1 opacity expansion
The zeroth order in opacity corresponds to the absence of medium effects, n(ξ) = 0, when no gluons are radiated, see equation (2.7). The first non-vanishing term in an opacity expansion of (2.2) is then the first order to absorb a factor (2π) 2 that is common to many formulas in the following.
In general, to any order in the opacity expansion of (2.2), factors |A(q)| 2 in the integrand appear always in the combination |A(q)| 2 − V totδ (q) . The terms V tot arise as a consequence of probability conservation, as we explain in more detail below.
To first order in opacity, see (2.16), the term proportional to V totδ (q 1 ) vanishes, and the result is of the form (2.10) of an elastic cross section |A(q)| 2 times a Bertsch-Gunion term (2.11). Hence, the N = 1 opacity contribution to the gluon energy distribution (2.2) accounts for all radiated gluons, which have interacted with exactly one scattering center in the medium. The prefactor (n 0 L) in (2.16) counts the number of independent gluon productions which occur within the length L.

N = 2 opacity expansion and formation time
Medium-induced quantum interference arises, if a single gluon is produced in interactions with at least two scattering centers. In the opacity expansion, this is realized for N ≥ 2. In particular, for N = 2, the medium-induced gluon distribution can be written in the form Here, we have adopted the following conventions [36]: To N -th order in opacity, subscripts are labeled such that i = 1 is the last, i = 2 the next to last and i = N the first scattering center along the trajectory of the partonic projectile. Also, the sign of the transverse momenta q i are chosen such that they are flowing from the projectile to the medium. The qualitatively novel feature of the N = 2 result (2.18), compared to the first order result (2.15), is the appearance of an interference factor In general, interference factors depend on the in-medium path length L and on transverse energies For the following, it will be useful to view the inverse of these transverse energies as formation times. In particular, The interference factor (2.19) interpolates between the two limiting cases 0 , for L τ 1 , n 0 L = const. (2.23) In both limiting cases, the energy distribution (2.18) has a probabilistic interpretation: • Incoherent production limit L τ 1 , n 0 L = const.
Here the Bertsch-Gunion term R(k + q 1 ; q 2 ) denotes a medium-induced radiation term, for which the gluon was produced incoherently on the first scattering center with momentum transfer q 2 and scattered incoherently on the last scattering center with momentum transfer q 1 .
Here, the Bertsch-Gunion term R(k; q 1 + q 2 ) denotes a coherent gluon production in which the two scattering centers are not resolved but act effectively as a single one.
In the expressions above, there are terms proportional to |A(q 1 )| 2 |A(q 2 )| 2 . These correspond to processes, in which the radiated gluon exchanges momentum with exactly two scattering centers. In addition, there are terms proportional to V tot |A(q 2 )| 2 , in equations (2.24) and (2.25), which involve only one momentum transfer with the target. For these latter terms, the totally coherent and incoherent limits differ by a factor 2. This can be understood in terms of a probabilistic picture of the partonic dynamics: In the incoherent case, the gluon can scatter on the second scattering center at ξ 1 only after it was produced incoherently at position ξ 2 . The corresponding weight from the integrals along the trajectory is ∝ L 0 dξ 2 n 0 L ξ 2 dξ 1 n 0 = (n 0 L) 2 /2. In contrast, in the coherent case when both scattering centers lie within the formation time of the gluon, their time ordering does not matter and the probability conserving contribution has the weight L 0 dξ 2 n 0 L 0 dξ 1 n 0 = (n 0 L) 2 , which is a factor 2 larger.

2.3
Combining medium-induced and vacuum gluon radiation (case ξ 0 = 0) In section 2.2, we discussed how destructive interference gives rise to formation time scales in the gluon energy distribution (2.2) with initialization ξ 0 = −∞, where vacuum radiation is absent. Here, we parallel this discussion for the initialization ξ 0 = 0, when the hard projectile splits also in the absence of a medium, as expected for a virtual state. To zeroth order in opacity, the gluon energy distribution (2.2) yields the singular part dI dz dk = α s C R 1 k 2 1 z . This is the leading order quark or gluon splitting function for z = ω/E within the eikonal approximation (2.1).
Destructive interference arises already to first order in opacity, (2.26) Here, τ 1 = 1/Q 1 is the formation time of the gluon prior to scattering on the medium with momentum transfer q 1 . The limiting cases are: • The limit L τ 1 , n 0 L = const. One finds the limit As we discuss in more detail in appendix A, this limit is consistent with the probabilistic picture that a gluon can only be produced in a scattering if it is formed as part of the incoming projectile wave function prior to the scattering.
On the right hand side of this equation, the first term is proportional to V tot and implements probability conservation: the total probability that a scattering with some momentum transfer q 1 occurs is subtracted from the N = 0-contribution that no momentum transfer occurs. If a momentum transfer occurs, then this momentum transfer can either shift probabilistically the transverse momentum of the fully formed gluon. This is the second term proportional to H(k+ q 1 ). Alternatively, the momentum transfer leads to a medium-induced gluon production, distributed according to the Bertsch-Gunion term R(k, q 1 ).

N = 2
Opacity expansion for ξ 0 = 0 For ξ = 0, the 2nd order in opacity of equation (2.2) can be written in the compact form [36] ω where (2.32) We consider again the case of a fixed number of effective scattering centers, n 0 L = const. In the limit n 0 L = const, L → 0, expression (2.29) vanishes, and so do all higher orders in opacity. In the opposite limit, n 0 L = const, L → ∞, one finds the totally incoherent limit The probabilistic interpretation of this expression is as follows: If the gluon has interacted incoherently with two scattering centers prior to escaping from the medium after length L with momentum k, then this gluon was either produced in a vacuum splitting and accumulated transverse momentum incoherently in two scattering. This is the term H(k + q 1 + q 2 ). Alternatively, the gluon was produced in a medium-induced interaction R(k + q 1 , q 2 ) with momentum transfer q 2 and accumulated additional transverse momentum q 1 incoherently in a second interaction. The second and third line of (2.34) readjust the probabilities that the gluon was produced with less than two momentum transfers from the medium. In particular, to all orders in N , the vacuum emission H(k) remains unmodified by the medium with the weight given by the no-scattering probability S = exp [−n 0 L V tot ], and the last line is the second order in opacity of S H(k). Similarly, the second line readjusts the probability for gluon production processes with exactly one scattering center involved. What dictates the scale at which the vanishing (totally coherent) radiation pattern (2.33) evolves into a fully developed incoherent radiation pattern (2.34)? For reasons that will become clear in the following subsection, we focus our discussion of this question on the medium-induced radiation term R(k + q 1 , q 2 ). We observe that in the limit n 0 L = const, L → ∞ of the gluon energy distribution (2.29), only the term proportional to Z 2 contributes to the medium-induced radiation term R. The limiting cases of Z 2 are lim n 0 L=const , L→∞ Inspection of equation (2.31) shows that for n 0 L = const, the first term vanishes for scales L 1/Q 2 and the second term for length scales L 1/Q 1 . To fully explore the physical implications of this observation, we recall that Q 2 is the transverse energy of the gluon prior to interacting with the target, and Q 1 is the transverse energy of the gluon after the first and prior to the second scattering. For the most likely scattering histories, transverse energy will be built up step by step in multiple scattering, Q 2 Q 1 . We have written this as a strong inequality with the idea that mediuminduced transverse momentum broadening should dominate over the initial transverse momentum of the vacuum radiation. Now, for Q 2 Q 1 , one sees that the second term in (2.31) dominates the value of Z 2 for sufficiently large L, and this second term dies out on length scales L 1/Q 1 . This leads us to the qualitative conclusion that it is the formation time 1/Q 1 of the gluon prior to its last interaction with the target that determines whether the radiation R takes place. The gluon is only radiated if its formation time is sufficiently short so that formation is completed on a scale comparable with the in-medium path length.

Guidance for an MC implementation
A remarkable simplification of MC simulations of the k-integrated radiation pattern arises from the fact that vacuum terms like H(k + q) in (2.28) do not contribute to parton energy loss. This is so, since H(k + q) amounts to a probability-conserving redistribution of gluons in transverse momentum space; this redistribution affects neither the yield of emitted gluons, nor their energy distribution. As a consequence, neglecting the terms proportional to H does not affect the gluon energy distribution ωdI/dω. For k-differential distributions, a similar a priori argument does not exist. We note as an aside that terms proportional to H were not taken into account in the original derivation of the BDMPS-Z formalism. They appeared first in the derivation of Ref. [36] that leads to (2.2). That they modify the transverse momentum distribution was also recognized in Ref. [46]. However, there is numerical evidence that inclusion of these terms is a numerically small effect [36]. Based on this observation, we shall seek a MC implementation of the BDMPS-Z formalism that neglects terms proportional to H. This treatment is exact for k-integrated quantities, and -as we shall show in section 6 -it is a satisfactory approximation for k-differential information.
For the medium-induced radiation terms R, at first order in opacity, the only difference between the cases ξ 0 = −∞ (2.16) and ξ 0 = 0 (2.26) is the reduction in the phase space of R due to the destructive interference term n 0 (L Q 1 − sin (L Q 1 )) /Q 1 . The analysis to first order in opacity did not allow us to disentangle between an interpretation of this phase space cut in terms of either i) the formation time prior to the very first or ii) prior to the very last interaction with the medium. The analysis of the 2nd order in opacity, however, gave support to the second interpretation, see section 2.3.2. Motivated by this observation, we shall propose in section 4 a Monte Carlo implementation of the BDMPS-Z formalism for ξ 0 = 0, according to which gluons are rejected from the simulation if their formation is not completed within the medium.
The analysis of the opacity expansion in section 2.3.2 supports only the parametric statement that those medium-induced gluons contribute to the distribution (2.2) whose formation is completed on a length scale comparable to L. It is one conceivable (though not unique) implementation of this parametric argument to count solely gluons whose formation is completed within the medium. We note that in establishing a one-to-one correspondence between the opacity expansion of (2.2) and a MC algorithm, this is the only point where we have found only parametric and not quantitative guidance. Accordingly, we have tested numerically some variations of this prescription, and we shall comment on this in section 5.

A simplified problem: a MC algorithm for N g in the totally coherent and incoherent BDMPS limits
The main aim of this paper is to formulate a MC algorithm that interpolates correctly between the analytically known BDMPS results in the opacity expansion. Explicit expressions for these limits are known analytically [36] to arbitrary high orders in opacity. For the case of an incident projectile (ξ 0 = −∞), the totally coherent limit is and the incoherent limit is Here, we have used the shorthand In general, contributions to N -th order in opacity contain products of a number N s (1 ≤ N s ≤ N ) of cross sections |A(q i )| 2 , and a number N − N s of cross sections V tot , obtained from expanding the prefactor exp [−n 0 L V tot ].
In this section, we consider first the simpler problem of formulating for the limits of totally coherent and incoherent gluon production an algorithm for the momentum-space integrated average number of radiated gluons, This study will be extended to the differential spectrum in section 4.

Relating BDMPS-Z to elastic and inelastic mean free paths
We consider first the N s = 1 scattering contribution to the totally coherent and incoherent BDMPS limits (3.1) and (3.2). The resulting average number of radiated gluons is Here, we have used the analysis of equation (2.16) to define the inelastic cross section for incoherent gluon production on a single scattering center as Here, the integrations over k and ω require regularization. The value of the regulator is a physical choice: it determines up to which soft scale infrared and collinear production processes are counted towards the inelastic cross section. We shall explain in section 5 how, based on this definition of σ inel , one can calculate measurable quantities that are insensitive to the choice of regulators. In the BDMPS-Z formalism, factors |A(q)| 2 and V tot are always multiplied by the density n 0 of scattering centers. The product n 0 σ inel defines the inelastic mean free path λ inel Physical results depend on λ inel , but they do not depend separately on σ inel and n 0 . As seen in the discussion of (2.3), the term |A(q)| 2 can be viewed as the differential elastic cross section dσ el dq for scattering of the partonic projectile on a single target. Accordingly, we identify The exponential factor exp (−n 0 L V tot ) can then be written in terms of the elastic mean free path λ el ,

Incoherent limit
The higher order terms of the coherent and totally incoherent BDMPS-Z limits (3.1) and (3.2) differ. In particular, for N s = 2, we have Here, the first term ∝ R (k, q 1 ) has a q 2 -independent integrand and can be written as a factor 1/λ el . This is a consequence of q 2 1 = V tot and the argument leading to (3.9). For the second term ∝ R (k + q 1 , q 2 ), a formal shift k → k − q 1 in the integral of (3.10) indicates that its contribution to the transverse momentum integrated average (3.10) is of the same magnitude. This prompts us to identify in the incoherent limit the higher orders of N s with Summing over all orders of N s , one finds This is the expected result for the average number of gluons produced incoherently within a length L, and it thus supports our identification of momentum-integrated terms in the BDMPS-Z formalism with elastic and inelastic mean free paths.

Totally coherent limit
To arbitrary order in opacity, we find from (3.1) for the totally coherent limit In general, the k-integration over R k, Ns j=1 q j differs from the integration over R k + j−1 l=1 q l , q j in the incoherent limit. However, both k-integrals are dominated by contributions from the two (IR regulated) singularities in the Bertsch-Gunion factor, and these dominant contributions are identical for both integrals. This prompts us to write (3.14) All totally coherent contributions N coh g (N s ) are exactly one factor 1/N s smaller than those of the incoherent limit (3.11). The resulting average number of gluons produced totally coherently is

Ambiguities in identifying mean free paths in the BDMPS-Z formalism
In the BDMPS-Z formalism, one calculates radiation cross sections for multiple scattering processes that have one additional gluon in the final state and that involve a very large number of elastic interactions. Therefore, the BDMPS-Z formalism is derived under the assumption that λ el λ inel .
The ratio of these mean free paths sets the value of the strong coupling constant, λ el /λ inel ∝ α s , see section 5 for a quantitative discussion. In this sense, the BDMPS-Z formalism is a weak coupling approach with regards to gluon radiation, whereas it resums the possibly non-perturbatively strong interactions between projectile and target. In general, the total mean free path λ tot is defined as However, in a formalism where λ el /λ inel = O(α s ) 1, the inverse of λ tot equals the inverse of λ el up to subleading corrections of O(α s ) that become negligible. That leads to some ambiguities in identifying mean free paths in the BDMPS-Z formalism. In the discussion so far, we have chosen to interpret V tot as a phase-space integrated elastic cross section. This is natural in the light of equation (2.3). On the other hand, one has also the choice of identifying n 0 V tot with 1/λ tot , and this ambiguity cannot be resolved within the accuracy of the BDMPS-Z formalism. We note that taking this alternative choice, one would find for instance N coh g = λtot λ inel 1 − e −L/λtot . In contrast to (3.15), this is smaller than unity for arbitrary values of λ inel and λ el , while equation (3.15) can be larger than unity for λ el > λ inel . In the region λ el λ inel , for which the BDMPS-Z formalism was derived, this difference becomes negligible.

MC algorithms for the incoherent and totally coherent BDMPS-Z limits
We consider a medium composed of scattering centers of a given density n 0 that provide elastic and inelastic cross sections to a projectile parton. We work within the approximations of the BDMPS-Z formalism, that means: We neglect elastic scatterings of the projectile partons, since they are unimportant for gluon radiation. And we neglect subsequent inelastic scatterings of the radiated gluons, since they are unimportant for understanding the energy loss of the projectile parton.

MC algorithm for the incoherent BDMPS-Z limit
We first formulate a MC algorithm that implements the BDMPS-Z formalism in the absence of quantum interference effects (incoherent limit). The starting point of the probabilistic evolution is a partonic projectile that propagates on a straight line ξ ∈ [0; L] through a medium of path-length L. The interaction between projectile and medium is characterized fully in terms of the inelastic mean free path λ inel of the projectile and the elastic mean free path λ el of the radiated gluons. The dynamic evolution starts at ξ = 0 and it proceeds according to the following steps: 1. Determine whether and where the projectile undergoes its next inelastic scattering Decide with probability 1 − S proj no (L) that a scattering occurs within the remaining in-medium path length L. Here, S proj no (L) is the probability that the projectile does not undergo any inelastic interaction within length L, If no further inelastic interaction is found, then stop the dynamical evolution. Else, determine the distance ξ to the next inelastic scattering center according to the probability density 2. After inelastic scattering, continue propagating the projectile After an inelastic interaction at position ξ, the outgoing projectile has a remaining in-medium path length L − ξ. To establish whether the projectile undergoes further inelastic interactions, repeat step 1 with inelastic no-scattering probability S proj no (L − ξ). Reiterate this step till no further inelastic interaction is found.

After inelastic scattering, propagate the produced gluon
The gluon, produced in an inelastic process at position ξ, has a remaining in-medium pathlength L − ξ. Determine the number and positions of additional elastic interactions of the gluon with the medium as follows: Determine whether and where the gluon undergoes its next elastic scattering, based on the elastic no-scattering probability S el That means, decide with probability 1 − S el no (L − ξ) that there is another elastic scattering, and determine its distance ξ − ξ according to the probability density Reiterate this process for each gluon till no further elastic scattering center is found.
According to this MC algorithm, the probability P inel (m) for generating dynamical scattering histories with exactly m inelastic interactions is determined by reiterating step 2 in the above algorithm, Since the algorithm produces exactly one gluon per inelastic interaction, P inel (m) is the probability for finding scattering histories with exactly m produced gluons. The average number of gluons per scattering history is which is consistent with the corresponding incoherent limit in the BDMPS-Z formalism, see (3.11).

MC algorithm in the presence of coherence effects
Coherence effects in gluon production processes can be accounted for by modeling the production as taking place over a finite formation time τ f in (2.21). The incoherent limit of gluon production is then realized for the case τ f λ inel , λ el and the totally coherent limit is realized for τ f L. To decide which of these limits applies to a specific gluon production process, the MC algorithm needs to know τ f . The dynamical determination of τ f requires k-differential information and will be discussed in the context of the k-differential algorithm in section 4. As a preparatory step, we explore here the formal limits τ f → 0 (incoherent) and τ f → ∞ (totally incoherent) gluon radiation, and we study in these limits k-integrated yields. We want to devise an algorithm that extends naturally to a k-differential version. To this end, we should use information about whether we work in the totally coherent or incoherent limit only in algorithmic steps in which information about τ f would be dynamically available in the k-differential version. Therefore, as long as the inelastic scattering and its kinematics is not yet determined, the MC algorithm must still allow for the cases that the inelastic production process could be either incoherent or could include coherence effects. This consideration prompts us to seek an MC implementation that starts from selecting an inelastic process as in the incoherent case, based on equations (3.17) and (3.18). Coherence effects will then be included by modifying the subsequent evolution and by reweighting the inelastic process that was selected with the probability of an incoherent production. Such reweighting is a standard Monte Carlo technique in algorithms that overestimate probabilities. We discuss now both these elements in more detail: Modifying the subsequent evolution: Assume that the MC algorithm has selected an inelastic process 'at ξ' according to (3.17), (3.18), and that the formation time τ f of the produced gluon is then found to be finite. How should this be taken into account in the further probabilisitic evolution? The general idea is that if τ f cannot be neglected ( τ f > λ inel , λ el ), then the position ξ selected in (3.18) cannot be interpreted as the 'point' of the gluon emission. Rather, we view the simulated pair of values ξ, τ f as specifying a region of extent τ f around ξ, over which the gluon production process takes place. Technically, this translates into the requirement that if gluon production could have started as early as ξ init = max [ξ − τ f ; 0], then the produced gluon is allowed to scatter elastically from time ξ init onwards, and not only after time ξ. This is a modification of step 3 of the incoherent algorithm. Physically, it means that within this entire region between ξ init and ξ init + τ f , elastic interactions act coherently with the inelastic one.
In the present subsection, we restrict our discussion to the totally coherent case, τ f L. In this particular limit, irrespective of the position ξ at which the MC algorithm allocates the center of an inelastic process, this process is delocalized over the entire in-medium path length L. As a consequence, irrespective of the choice of ξ, the radiated gluon can accumulate additional elastic interactions between ξ init = 0 and L.

Reweighting inelastic processes:
In the incoherent case, the probability that the projectile parton undergoes one or more inelastic interactions is given by 1 − S proj no (L), see (3.17). Each scattering center serves as an independent source of gluon production. In contrast, in the presence of coherence effects, it is the ensemble of several scattering centers that acts effectively as one source of gluon production. Therefore, the factor 1 − S proj no (L) overestimates the probability of inelastic interactions, and a reweighting is needed.
To determine this reweighting factor, we observe that in the totally coherent limit of the BDMPS-Z formalism, N coh g (N s ) in (3.14) denotes the average number of gluons produced with exactly (N s − 1) elastic and one inelastic interaction. The corresponding expression in the incoherent limit is given in (3.11) and it is one power of N s larger, N incoh g (N s ) = N s N coh g (N s ). Therefore, the N s -averaged number of emitted gluons can be obtained in the totally coherent limit, if a gluon selected according to (3.17), (3.18) and having undergone N s scatterings is accepted with probability Based on these considerations, we propose the following MC algorithm for the totally coherent limit: 1. Determine whether the projectile undergoes an inelastic scattering.
As in the incoherent case, use (3.17) to decide with probability 1 − S proj no (L) that a scattering occurs within the in-medium path length L. If no inelastic interaction is found, then stop the dynamical evolution.
Establish whether the projectile undergoes further inelastic interactions by searching with probability 1 − S proj no (L − ξ) for further inelastic scatterings between ξ and L.
3. After inelastic scattering, propagate the produced gluon up to length L and reweight its production probability.
In the totally coherent case, the production is delocalized over the entire medium of length L and therefore, all gluons undergo elastic scattering over an in-medium path length L. With probability 1 − w = 1 − 1 Ns , the produced gluons are rejected.

Validating the proposed MC algorithms
We have written MC programs that implement the algorithms proposed in sections 3.2.1 and 3.2.2 for the case of incoherent and totally coherent gluon production, respectively. To check that these algorithms reproduce the analytically known results of the BDMPS-Z formalism, we establish that they account for the average number of gluons produced per scattering history in both limits, N coh g and N incoh g . In addition, the MC algorithms allow us to plot the average number of gluons N g (incoh) j and N g (coh) j , produced with exactly j momentum transfers from the medium. Here, we test against this more differential information.
In the totally coherent limit, we see from equation (3.13) that the expansion of N coh g (N s ) to order N s involves gluon radiation terms with exactly N s momentum transfers. As a consequence, the average number of gluons produced with exactly j momentum transfers is given by (3.23) The analogous identification of orders in the opacity expansion with number of momentum transfers does not hold in the incoherent limit. As one sees for instance from equation (3.10), the second order receives contributions from gluons that were produced either with one single momentum transfer (these are the terms R(k, q 1 )) or with two momentum transfers (these are the terms R(k + q 1 , q 2 )).
To identify all contributions with a fixed number of momentum transfers, we write the incoherent limit of the BDMPS-Z formalism as a series Here, contributions involving the radiation term R k + j i=2 q i , q 1 denote gluon production processes with j-fold scattering (i.e. with (j − 1)-fold elastic scattering). Integrating formally over phase space, one finds that the average number of such gluons per event, produced with j-fold scattering, can be expressed in terms of complete and incomplete Γ-functions, One can check that the average number of incoherently produced gluons is again given by

A k-and ω-differential MC algorithm in the totally coherent and incoherent BDMPS-limits
In the previous section, we have shown how the coherent and incoherent limits of the phase space integrated average number of radiated gluons dω dk dI dω dk can be simulated in a probabilistic MC algorithm. In the present section, we extend these algorithms to a simulation of the differential gluon distribution dI dω dk . The basic building block for the differential distribution dω dk dI dω dk is the inelastic interaction of the projectile with a single scattering center. According to eqs. (2.10) and (2.16), the corresponding inelastic cross section is (4.1) We seek a MC algorithm that interpolates between the coherent and incoherent limits by treating all momentum transfers during the formation time of a gluon as coherent, and all scatterings outside the formation time as incoherent. Such an algorithm must keep track of the kinematics of the scatterings, and it must account dynamically for changes in the formation time. We propose an algorithm that as criterion for decoherence of the gluon requires the relative phase of the radiated gluon to become unity. More precisely, we observe that the interference factor (2.23) extracted from the BDMPS-Z formalism is best approximated by a Θ-function of the form 1 Therefore, we define the formation time by the condition We first discuss in section 4.1 the inputs and approximations of (4.1) that simplify an MC implementation. We then specify a MC algorithm before discussing how some of these approximations can be relaxed.

Inputs and approximations in the proposed MC algorithm
In the study of parton energy loss models and the BDMPS-Z formalism, a standard parametrization of elastic scattering cross sections is in terms of a Yukawa potential with a screening mass µ, (4.5) In the following, we work with this ansatz for µ ∈ [100 MeV; 1 GeV]. In equation (4.1), the inelastic cross section for a single incoherent scattering factorizes into the product of the elastic cross section and a radiation term. The term R(k; q) specifies how gluons produced with energy ω are distributed in transverse phase space prior to undergoing subsequent interactions. What matters for the decoherence of the gluon is its relative momentum with respect to the outgoing parent parton. If the final transverse momentum of the gluon is build up by many interactions with the medium, then the precise distribution of the transverse momentum at the inelastic interaction can be expected to be unimportant. Moreover, even if there are not many interactions with the medium, the transverse momentum of the gluon at the inelastic interaction 1 The interference factor f (x = L Q 1 ) = 2 (1 − cos x) /x 2 decreases continuously from f (0) = 1 to f (2π) = 0, and it oscillates for larger values of x with rapidly decreasing amplitude ∝ 1/x 2 . One finds will be set by the recoil received by the medium. These considerations prompt us to adopt the following approximation that simplifies the numerical implementation significantly In section 5, we shall provide numerical evidence that the approximation (4.6) is sufficient for a quantitative MC implementation of the BDMPS-Z formalism. With the help of (4.6), the total inelastic cross section simplifies to Here, we have considered gluon radiation in the range ω ∈ [ω min ; ω max ]. We note that the first line of (4.6) needs to be regularized, since the integral over R(k; q) is infrared divergent. Performing the integral over R(k; q) with an infrared cut-off around k = 0 and k = q, one finds f prop = 2 π [log (µ 2 / ) + const.]. In our MC algorithm, the infra-red regulator will not appear. Rather, for one arbitrary choice of model parameters, we shall adjust f prop such that the BDMPS result for the average parton energy loss is reproduced with the correct norm. For all other choices of model parameters, f prop is then kept fixed and the MC algorithm generates normalized results. What can be said a priori about the numerical value of f prop is that there is no physical reason for choosing an infrared regulator that is much smaller than the momentum scale µ. Therefore, the logarithm log (µ 2 / ) should not be large, and f prop should be of order unity. We shall confirm this expectation in section 5. We pause to comment on this approximation from a wider perspective: The BDMPS-Z formalism (2.2) does not depend on total elastic and inelastic cross sections, but only on the dipole cross section (2.4) that does not require regularization since it is differential in configuration space. However, the opacity expansion of (2.2) rearranges this formalism in a series that does contain total phase-space integrated quantities. To arrive at a probabilistic implementation, we have assigned to some terms in the opacity expansion of (2.2) the natural physical meaning of elastic and inelastic cross sections and of mean free paths (see eqs. (3.7) and (3.9)). This can only be done with the help of approximations and regularizations that are not explicit in the BDMPS-Z formalism (2.2). For instance, the identification of phase-space integrated expressions of the opacity expansion with rational functions of mean free paths (such as e.g. eq. (3.10)) is strictly speaking a physically motivated assignment rather than an analytically derived fact, since the transverse momentum integrals are infra-red divergent. The crucial test for the MC implementation is then that physical results do not depend on the regularization prescriptions employed and that they account quantitatively for the BDMPS-Z formalism (2.2). That this is so will be demonstrated in section 5.

A k-and ω-differential MC algorithm interpolating between the incoherent and totally coherent BDMPS-limits
1. Initialisation Set remaining path length of the projectile to total path length, L proj = L.

Determine whether and where the projectile undergoes its next inelastic scattering This step is implemented as described by equations (3.17), (3.18) and accompanying text.
If an inelastic scattering is generated at position ξ, then the remaining path length of the projectile is set to L − ξ. The produced gluon is propagated further according to the step 3 below. The algorithm repeats step 2 till no further inelastic scatterings are found in the remaining path length.

Kinematics of gluon emission and dynamical evolution of formation time
In the BDMPS-Z formalism, the gluon energy is distributed according to 1/ω. From this distribution, the gluon energy is generated. The initial transverse momentum of the gluon is generated from the distribution |A(k)| 2 ; the initial gluon phase is taken to vanish, ϕ = 0; the number of momentum transfers to the gluon is set to N s = 1, and the initial formation time is determined according to Then set the remaining gluon path length to the total path length, L gluon = L, and check for further elastic momentum transfers within the formation time: • With probability 1 − S el no (min(τ f , L gluon )) there is one more scattering. Determine the distance ∆L to the scattering centre and update the path length, L gluon = L gluon − ∆L, and the gluon phase Determine the momentum transfer q Ns from the scattering centre according to |A(q Ns )| 2 , set the transverse momentum of the gluon to k = Ns i=1 q i , and set N s = N s + 1. Iterate this point until no further scattering is found.
• With probability S el no (min(τ f , L gluon )) there is no further scattering. Continue with point 4.
4. Reweight the gluon production probability, and propagate gluons further. The gluons simulated in point 3 are trial gluons that have been selected with an overestimated production probability. Reweighting is needed to correct for this overestimate. If a trial gluon is generated with N s scattering centers within its formation time, then • With probability 1 − 1/N s , reject the gluon from the sample.
• With probability 1/N s , accept the gluon as part of the scattering history. Determine the end of the formation process of the gluon by localizing a formation time interval τ f in an arbitrary fashion around the initial production point ξ. Then determine further elastic momentum transfers to the gluon within the in-medium path length after formation has been completed. (This last step is needed only for the simulation of k-differential spectra.)

Accept only medium-induced gluons
To reproduce the radiation spectrum (2.2) for ξ 0 = 0, accept only gluons that are fully formed prior to leaving the medium.
It is a consequence of the approximation (4.6), that the gluon transverse momentum is build up identically in the coherent and incoherent case. We note as an aside, that it is possible to amend the above proposal such that it does not invoke the approximation (4.6). To do so, one has to start from the observation that a gluon produced with N s coherently acting scattering centers is produced according to the probability In our simplified algorithm, this expression is approximated by a factor λ el /λ inel , and λ inel specifies the probability with which an inelastic scattering occurs. There are standard reweighting techniques that would allow one to overestimate the probability of inelastic interaction and to then correct it to the factor (4.10). In the present work, we did not exploit this numerically more demanding procedure, and we did not find any indication that such a procedure is needed to reproduce quantitatively the BDMPS-Z formalism (2.2). The idea that the concept of formation time plays a central role in the probabilistic implementation of medium-induced gluon radiation has been formulated previously. However, in our effort to arrive at a quantitatively reliable, probabilistic, formation time based formulation of the BDMPS-Z formalism, we had to overcome several conceptions that were naively assumed at least by us, but possibly also by others. In particular, a MC formulation that selects gluon production processes according to an incoherent inelastic scattering probability overestimates gluon production in the presence of interference effects. A quantitatively reliable implementation must correct for this overestimate, and the algorithm proposed here is, as far as we know, the first one that does so. On general grounds, one expects that this feature is not specific for the BDMPS-Z formalism, but persists in more complete formulations of radiative parton energy loss. Secondly, it turns out that the BDMPS-Z formalism cannot be implemented exactly in a formulation that interprets formation times as deadtimes for subsequent gluon production. Technically, this can be seen from the form of the average number of radiated quanta N g j as a function of the number of active scattering centers j, discussed in subsection 3.3. (Formulations based on a dead time interpretation would lead to expressions for N g j that contain terms ∝ λ inel in the arguments of exponentials.) That formation times are not dead times for subsequent gluon production could have been expected on the simple ground that the BDMPS-Z formalism is based on a multiple scattering calculation with only one gluon in the final state and therefore cannot account for the destructive interference between different gluons. It remains to be seen whether this feature persists in more complete analytical calculations of medium-induced gluon emission.

Numerical results on the gluon energy distribution
The MC algorithm of section 3 and 4 is tailored to provide a probabilistic implementation of the opacity expansion of (2.2). At fixed order in opacity, terms in (2.2) can be pictured as arising from interactions of the partonic projectile with a fixed number of scattering centers. This discrete picture of the medium lends itself naturally to a MC implementation, and the proposed algorithm reproduces the analytically known distribution in the number of scattering centers, see Fig. 1.
In contrast, in the multiple soft scattering limit of (2.2), information about the discrete structure of the medium is lost. This limit is obtained from a saddle point approximation of the path integral in (2.2), setting n σ(r) = 1 2q r 2 . In this approximation, the BDMPS-Z transport coefficient q characterizes the average transverse momentum squared, transferred from the medium to the projectile per unit path length. The medium can be pictured as providing for the projectile a continuous transverse color field whose strength is characterized byq.
Here, we shall compare results of the proposed MC algorithm to the BDMPS-Z multiple soft scattering approximation of (2.2) according to which the energy distribution (2.2) of gluons emitted from a highly energetic projectile shows the characteristic 1/ √ ω-dependence of the non-Abelian Landau-Pomeranchuk-Migdal effect, This 1/ √ ω-spectrum is cut-off due to formation time effects at a characteristic gluon energy ω c = 1 2q L 2 . Integrating ω dI dω , one finds the average parton energy loss Here, the critical path length L c is the maximal coherence length, which occurs for the maximal kinematically allowed gluon energy ω max (typically taken to be the projectile energy E proj ). For lengths L > L c , one expects hence that different regions of the medium act incoherently to gluon production and that ∆E(L) increases linearly with L. The differential distribution (5.1) continues to show the characteristic coherence effects for L > L c , since each gluon entering this distribution was produced coherently over a distance τ f that depends on ω.

Multiple soft scattering limit in the MC algorithm
To realize the multiple soft scattering approximation in the MC algorithm, we ensure first that there are many elastic interactions per inelastic mean free path. Hence, we shall work in the limit Moreover, we ensure that all elastic interactions are soft by cutting off the power-law tails of the Yukawa scattering potential (4.4) at |q| = 2 µ, This approximation in the MC algorithm can be shown to correspond on the analytical side to a saddle point approximation of the path integral (2.2) by writing in (2.4) σ(r) ∝q r 2 .
The soft multiple scattering approximation of (2.2) and the average parton energy loss (5.2) are functions of α s C R and forq, which are not input parameters of the MC simulation. Instead, one specifies for MC simulations the elastic and inelastic mean free paths, and the average transverse momentum transfer µ in the elastic scattering cross section. To express the BDMPS energy loss formula in terms of these input parameters, we rewrite the strong coupling constant with the help of eq. (4.7), .
From the MC simulation, we determine the event averaged squared transverse momentum q 2 transfered from the medium to a radiated gluon per unit path length L p , We then define operationally 2q In general, µ 2 /λ el would be a poor approximation of q eff , but for the particular choice of soft scattering centers (5.4) regulated at |q| = 2 µ, q 2 = µ 2 and q eff agree with µ 2 /λ el . We can now express the BDMPS parton energy loss formula (5.2) in terms of input parameters of the proposed MC algorithm, It is this form of the BDMPS parton energy loss formula that we test in the MC studies presented in this section. In the following subsections 5.2 and 5.3, we explore the proposed MC algorithm for values λ el O(10 −1 ) λ inel that realize the multiple scattering approximation (5.3). We note that the strong coupling constant in (5.3) is proportional to λ el /λ inel ; moreover, it decreases with a large logarithm 1/ log (ω max /ω min ) O(10 −1 ). (Unless stated otherwise, the numerical results in this section are for ω max = 100 GeV and ω min = 50 MeV.) As a consequence, the numerical values for the average energy loss presented in the next subsections 5.2 and 5.3 will be typically a factor 10 −2 lower than realistic values, since they have been obtained with an artificially low strong coupling constant. It is only by relaxing the multiple soft scattering approximation (5.3) that realistic values of the strong coupling strength can be implemented in the present MC algorithm. This will be done in section 5.4.
2 In a simplified scenario in which a fixed k 2 = µ 2 is transferred per mean free path λ el from the medium to a gluon, the MC algorithm will accumulate within a length L = n λ el a gluon phase ϕ ≈ 1 2ω n−1 j=0 j µ 2 λ el 1 2ω This phase differs by a factor 2 from the standard analytical pocket estimate ϕ = k 2 2ω L 1 2ωq L 2 . The reason is that the squared transverse momentum k 2 ∆L accumulated between L − ∆L and L, can contribute to ϕ only with k 2 ∆L ∆L/2ω and not with k 2 ∆L L/2ω. This illustrates that pocket formulas for ϕ (and a fortiori for ω c and L c ) should not be expected to provide numerically accurate prefactors but identify the parametric dependencies only.  Fig. 2 shows the medium-induced gluon spectrum for a projectile parton propagating through a medium of path length L. These and the following results were obtained for MC simulations of N evt = 10 6 events. For sufficiently large in-medium path length L, the spectrum ω dI dω approaches the characteristic 1/ √ ω-dependence expected for the non-abelian LPM effect. This dependence may be understood also by the following parametric argument: In the incoherent limit, gluon production on a single scattering center results in a spectrum ∝ 1/ω. Coherence effects imply that the number N coh of scattering centers located within the formation time of the gluon act as one single effective scattering center. The resulting gluon spectrum is ∝ 1 N coh ω . The average number of coherently acting scattering centers is proportional to the average formation time, and this average formation time should satisfy t coh ∝ ω q t coh . As a consequence, N coh ∝ t coh ∝ √ ω and therefore coherence effects change the gluon energy spectrum by one factor √ ω. For sufficiently small in-medium path length L or sufficiently large projectile energy ω max , the formation of gluons of high energy ω is suppressed since their formation time becomes comparable to the entire in-medium path length. Parametrically, this suppression is expected to set in at a characteristic gluon energy ω c = 1 2q L 2 , that takes values of ω c = 10, 40 and 90 GeV respectively for the in medium path lengths L = 1, 2 and 3 fm explored in Fig. 2. We note as an aside that in the limit ω c L → ∞, the expression (2.2) reduces to the BDMPS limiting result dI dω ∝ log cos 1+i √ 2 ωc ω . Numercial inspection of this limiting case reveals that the transition from the small-ω to the large-ω behavior of (5.1) occurs at values that are a factor ∼ 3 − 5 smaller than ω c . This is quantitatively consistent with the location of this transition region in Fig. 2, and it further illustrates the comment in footnote 2. Furthermore, a lower value for the transition energy was also found in [47,48]. We conclude that the proposed MC algorithm reproduces the ω −3/2 -dependence of the BDMPS-Z formalism for soft gluon production up to the expected scale which is of order ω c . For higher gluon energies, one observes a steeper ω-dependence, consistent with the BDMPS-Z formalism, but one finds some deviations from the ω −3 -dependence of (5.1) for realistic projectile energies. Since gluon energies ω ω c are known to be numerically unimportant in the BDMPS-Z formalism, these deviations will turn out to be negligible for the following.

MC results of the gluon energy distribution and control of cut-off dependence
We now turn to an issue that is crucial for the predictive power of a MC algorithm, namely that the physics results of the algorithm are insensitive to the numerical choices of IR and UV regulators, though various intermediate steps in the algorithm may depend on the choice of such regulators. To be specific, the MC algorithm selects inelastic interactions with a probability 1 − exp (−L/λ inel ) that depends on the total inelastic cross section. This cross section (4.7) depends explicitly on IRand UV regulators ω min and ω max . The physics output will still be insensitive to these regulators if the dependence of the total inelastic cross section on phase space available for radiation is respected in the MC implementation. Technically, this is achieved in the present algorithm by rescaling λ inel according to the cut-off dependence of σ inel . Fig. 3 illustrates that with this rescaling, the proposed MC algorithm satisfies this important property of cut-off independence. More explicitly, by changing the values of the IR and/or UV cut-off, we change the numerical value of σ inel so that λ inel varies like λ inel ∝ 1/ log (ω max /ω min ). Once an inelastic scattering center is identified in the MC simulation, the kinematics of the emitted gluon is then chosen in the same kinematic range ω ∈ [ω min ; ω max ] that was used for the calculation of σ inel . As seen on the right hand side of Fig. 3, this procedure results in cut-off independence of physical results: choosing ω min and ω max specifies the range within which results are generated, but it does not affect the results within this range.
In general, the appearance of IR and UV cut-offs in the MC algorithm can have different reasons. For differential inelastic cross sections that implement exact energy-momentum conservation, there is no need to specify by hand an UV cut-off ω max . Rather, the form of the cross section will automatically account for the physical requirement that gluons can only be emitted with energies smaller than the energy of the incoming partonic projectile, ω max = E proj . The introduction of an UV cut-off is only necessary, since one uses typically the approximate high-energy limit of the radiation cross section ∝ 1/ω, that extends to arbitrarily large gluon energy. For the case of the IR cut-off ω min of the ω-integration, or for the case of the corresponding IR regulator of the k- differential cross section that enters the total inelastic cross section(4.7) via the factor f prop , the situation is different. There is no perturbative physics argument that could specify the precise value of these regulators. All one can require is that whatever values for these IR cuts are chosen, the physics simulated above these values does not depend on the choice of the regulator. Fig. 3 illustrates that this requirement is satisfied by the proposed MC algorithm.

MC results for the average parton energy loss ∆E
In this subsection, we discuss MC simulation results for the average parton energy loss ∆E. The main purpose of this discussion is to give numerical support to equation (5.8). Fig. 4 shows the L-dependence of the average energy loss for different values of the UV regulator ω max = E proj . In this and all subsequent simulations, the value of the inelastic mean free path was adjusted to the varying phase space, λ inel ∝ 1/ log (ω max /ω min ), so that cut-off independent results were obtained. Since gluons of larger energy ω require on average a longer in-medium path length to form, one expects on general grounds that the small-L behavior of the average parton energy loss ∆E is independent of the choice of the UV regulator ω max . This is seen clearly in Fig. 4 for sufficiently small L. Moreover, for sufficiently large L > L c , Fig. 4 confirms the expected linear Ldependence of ∆E. The transition from a quadratic to a linear dependence occurs at an in-medium path length of order L c ∝ √ ω max that increase with the UV cut-off ω max = E proj . Our remark about the accuracy of scale estimates, made about ω c in the discussion of Fig. 2, and in footnote 2, applies here too. We note in particular that the quantity L c is not a quantitative prediction of the BDMPS-Z formalism, but that it characterizes only the expected parametric dependencies of (2.2). Consistent with this, we observe that the transition from quadratic to linear behavior shows the parametric dependencies expected for the BDMPS-Z formalism.
In general, we find that the ansatz ∆E(L) = a L 2 provides a very good description of results of the MC simulations, if the prefactor a is fit in the range L < L c . But for values L L c , results for ∆E(L) tend to lie significantly below the L 2 -fit. As we discuss now, this deviation can be understood by studying the dependence of average parton energy loss on the elastic mean free path λ el , see Fig. 5. According to equation (5.8), ∆E ∝ λ el q eff . Since q eff =q ≈ µ 2 λ el , one expects that the average parton energy loss is independent of λ el for L < L c and for fixed average momentum transfer µ per scattering center. On the other hand, the critical length depends on λ el , L c ∝ 1/ √q ∝ √ λ el and therefore the L 2 -dependence of ∆E should extend to larger values of L for larger values of λ el . On the scale of sufficiently large L, these features are confirmed by the MC data, see the right hand side of Fig. 5: results fall on a common L 2 -curve for L < L c , and they turn to a linear L-dependence at a scale L c ∝ √ λ el . (The curve for λ el = 0.01 fm in Fig. 5 turns to a linear L-dependence around L ∼ 5 fm, while the curve for λ el = 0.05 fm shows a quadratic behavior to much larger L, data not shown.) On scales of very small in-medium pathlength, however, there is a characteristic deviation from the λ el -independence of ∆E. On the left hand side of Fig. 5, we fit an L 2 -dependence to the data obtained for the smallest elastic mean free path λ el = 0.002 fm. Remarkably, at large L, this fit reproduces perfectly the data simulated with a 25 times larger mean free path, although this parameter set lies significantly below the L 2 -curve for L < 1.0 fm. This illustrates that increasing λ el at fixed L amounts to studying deviations from the soft multiple scattering limit. These occur when the probability for no scattering becomes sizeable -this is an effect that does not occur in the multiple soft scattering calculation but that will always be present in the Monte Carlo implementation. However, at fixed small value L, the characteristic L 2 -dependence of the soft multiple scattering limit (5.8) can always be recovered by going to sufficiently small values of λ el . The default parameter choice λ el = 0.01 fm used in this section was largely motivated by the idea to go sufficiently deep into the multiple scattering limit λ el λ inel to observe a quadratic L-dependence on a scale of 1 fm. In summary, Fig. 5 confirms the λ el -dependence of the expression (5.8) for the average parton energy loss and it quantifies the relevance of the multiple soft scattering approximation (5.3).
Motivated by these observations, we perform all fits of the quadratic L-dependence of ∆E(L) in the range L ∈ [0; L c ]. We can then confirm the other parametric dependencies of equation (5.8).
In particular, we have calculated the average parton energy loss for different values of the inelastic mean free path λ inel , and we have fit the prefactor a of ∆E(L) = a L 2 , see Fig. 6. The average energy loss is found to be inversely proportional to λ inel .
We have also characterized the dependence of the average parton energy loss on the quenching parameter. The right hand side of Fig. 7 confirms that for the current choice of soft elastic scatterings (5.4), the effective quenching parameter q eff satisfies indeed q eff =q ≈ µ 2 /λ el . The left hand side of Fig. 7 provides the check that the average parton energy loss depends linearly on q eff .
With the figures 4, 5, 6 and 7, we have confirmed all parametric dependencies of the BDMPS-Z result (5.8) for the average parton energy loss. To determine the overall normalization, we make the ansatz ∆E = c q eff L 2 . A direct fit of ∆E = c q eff L 2 to MC data in Fig. 7 results in c = 0.0021. On the other hand, we require from (5.8) For the values λ el /λ inel = 1/10, and log (ω max /ω min ) = log(100/0.05) ≈ 7.6 used in the simulations of Fig. 7, we find therefore f prop = 1.58 . (5.10) We recall that in the present formulation, the value of f prop is not a prediction of the BDMPS-Z formalism. Rather, as argued in the discussion of (4.7), this factor absorbs the remaining dependence on the IR cut-off of the total inelastic cross section that is needed in intermediate steps of the MC algorithm. It is a prediction, however, that the factor f prop is of order unity, and that it is a universal factor that is valid for all model parameter choices. This later statement will be further supported by the numerical studies in section 5.4. We note that the factor f prop absorbs also uncertainties of the MC implementation. In particular, we know from further numerical studies that f prop grows roughly proportional with ϕ(τ f ) (data not shown). Since the choice of ϕ(τ f ) = 3 adopted here is uncertain by ca. 15% (see discussion of eq. (4.3)), the factor f prop will also absorb this uncertainty. Noting that in the proposed MC algorithm, the acceptance criterion for produced gluons (step 5 in section 4.2) is based solely on the parametric arguments of section 2.3.3, we have also investigated modifications of this acceptance criterion. In one extreme alternative version, we required instead that gluons are counted towards the medium-induced spectrum if their formation after the last momentum transfer is completed within a time of scale L, irrespective of whether this amounts to completed formation inside or outside of the medium. For this modified MC algorithm, we repeated the entire study of sections 5 and 6 with analogous conclusions and very similar figures. The main difference compared to the results presented here was that we found an f prop that was approximately a factor 2 smaller than the value quoted in (5.10). From this we conclude that depending on how one implements those elements of the MC algorithm for which the opacity expansion of (2.2) provides only qualitative but not quantitative guidance, one arrives at a different factor f prop of order unity. Most importantly, however, once these ambiguities in the MC implementation are fixed by choosing a specific value for f prop , the absolute normalization of the simulated parton energy loss is fixed for all parameter choices.

MC results for phenomenologically motivated parameter values
The choice of elastic and inelastic mean free paths amounts to specifying the strong coupling constant α s , see (5.5). In the numerical studies in subsections 5.2 and 5.3, we focussed on the perturbative limit λ el λ inel . The parameters chosen in these studies correspond to a nominally perturbative regime in which α s ∼ O(10 −2 − 10 −3 ) or smaller. We now establish that the properties of the MC algorithm discussed in sections 5.2 and 5.3, persist for phenomenologically more relevant parameter choices. According to equation (5.5), realistic values for α s are obtained for choices λ el = O(λ inel ), and this motivates the parameter choices of the simulations shown in Fig. 8 and 9. These simulations included gluon radiation in the range ω min = 50 MeV to ω max = 100 GeV. We note that physical results do not depend on the precise choice of ω min ; in particular, a larger value of ω min could be absorbed in a rescaled inelastic mean free path λ inel ∝ 1/ log (ω max , ω min ), as discussed in the context of Fig. 3. On the other hand, physical results depend on the upper boundary ω max that sets the critical length L c 4ω max /q at which ∆E(L) changes from a quadratic to a linear L-dependence.
The scale of ω max is set by the physical UV cut-off on the radiation spectrum that is given by the energy of the partonic projectile. For the simulations shown in Fig. 8, we studied two different values of λ inel = λ el for Yukawa masses µ = 0.2, 0.5, 0.7 and 1.0 GeV in the elastic scattering cross section (4.5), respectively. These Yukawa masses set the scale of the transport coefficient q eff ∼ µ 2 /λ el . With these parameter choices, medium-induced gluon radiation is studied for a projectile parton of E proj = 100 GeV energy, propagating through a time-independent static medium of transport coefficient q eff . The results in Fig. 8 shows that the transition from a quadratic to a linear behavior occurs also for phenomenologically relevant parameter values at a scale of order L c , as established in section 5.2 in the multiple soft scattering limit. Fitting a quadratic dependence to the small-L region of ∆E(L), we confirm all parametric dependencies of (5.8). Moreover, we confirm within an accuracy of better than 5 %, that the proportionality factor f prop of (5.8) takes the same value as in the multiple soft scattering limit, f prop = 1.58. This shows that for one universal normalization of the inelastic cross section (4.7), the MC algorithm accounts faithfully for the results of the BDMPS-Z formalism (2.2) over a very wide parameter range, including phenomenologically motivated parameter choices.  The BDMPS-Z path integral (2.2) does not depend separately on the coupling constant, the number of scattering centers per unit path length n and the dipole cross section σ. Rather, it depends only on α s and on the linear combination n σ. As a consequence, the MC implementation of (2.2) does not depend separately on q eff , λ el and λ inel . Rather, it depends only on two combinations of these three parameters, which may be chosen to be q eff and λ el /λ inel , only. This is clearly seen in Fig. 8, where the choices λ el = λ inel = 0.1 fm with µ = 0.7 GeV and λ el = λ inel = 0.2 fm with µ = 1.0 GeV correspond to different microscopic pictures of the interaction between projectile and medium, but result both in the same average squared momentum transfer per unit path length q eff = 5 GeV 2 /fm, and in the same average parton energy loss.
Comparing Fig. 9 to Fig. 2, we observe that also the ω-differential information continues to show for phenomenologically motivated parameter choices the main features that we have established in the multiple soft scattering limit. In particular, the spectrum shows for soft gluon energies the ω −3/2 -dependence characteristic for medium-induced coherence, and for large gluon energy a steeper fall-off ∝ 1/ω 3 . Also, the transition between these two limiting spectra occurs at the scale ω ∼ ω c = 1 2q L 2 , as expected from the BDMPS-Z result (5.1). We finally note that for the model parameters chosen in this subsection, one finds the still rather small value of α s ≈ 0.1. For larger values of α s , the resulting average parton energy loss will increase correspondingly. Therefore, Fig. 8 illustrates that for phenomenologically relevant parameters and length scales, the average parton energy loss can attain values of tens of GeV.

Numerical results on transverse momentum broadening
In section 5, we have demonstrated that the MC algorithm of subsection 4.1 reproduces faithfully the ω-dependence of (2.2) for k-integrated quantities. We now discuss how the MC algorithm accounts for the k-dependence of the BDMPS-Z formalism.
It is a generic feature of the BDMPS-Z formalism that the transverse momentum of produced gluons is accumulated according to transverse Brownian motion, k 2 ∝q L . (6.1) To identify this feature numerically, we plot in Fig. 10 the simulated double differential distribution dI dω dk for different ranges of gluon energies ω as a function of κ 2 = k 2 /q eff L. In accordance with (6.1), the main contribution to the yields of simulated gluons lies in the range κ 2 ≤ 1, irrespective of the gluon energy, and irrespective of the choice of the parameters λ el , λ inel and µ 2 that control the rate of gluon production and its transverse momentum broadening. The double differential distribution ω dI dω dκ 2 of (2.2) has been analyzed and plotted for the soft multiple scattering limit and the N = 1 opacity approximation in Ref. [47]. We note that the results of the MC simulation shown in Fig. 10 reproduce very well the main results of Ref. [47]. In particular, the gluon yield dies out at a scale κ 2 ∼ O(1), it decreases with increasing gluon energy, and it shows a plateau for logarithmically small values of κ 2 . Also, the overall normalization of the MC results for the double differential distribution is in general agreement with the results of Ref. [47].
There are also qualitatively noteworthy though quantitatively small differences between the analysis of (2.2) in Ref. [47] and the output of the MC algorithm proposed here. In particular, destructive medium-induced interference effects can modify the gluon radiation such that in comparison to the vacuum distribution, the total yield of produced gluons is reduced in some phase space region below its average value in the vacuum. This would show up in negative values of the medium-induced gluon energy distribution ω dI dω dk . Such an effect has been identified indeed in a small phase space region of (2.2) [47]. A similar observation of small negative contributions has been made for the k-integrated distribution ω dI dω at very small in-medium path length. In contrast to the analytic calculation, which considers both vacuum and medium induced radiation and subtracts the unperturbed vacuum spectrum from the total gluon spectrum, the MC algorithm neglects the vacuum emissions (as a consequence the MC spectrum cannot become negative). While this can be seen in small deviations of (2.2) from MC results, we emphasize here that all numerically important, generic features of (2.2) are accounted for quantitatively by the MC algorithm.
It is also a characteristic feature of the BDMPS-Z formalism that medium-induced gluon radiation occurs for all gluon energies that accumulate at least a phase factor of order unity in the medium, k 2 L 2 ω = ωc ω > 1. In combination with the transverse momentum broadening (6.1), the technical manifestation of this statement is that the gluon radiation spectrum is only a function of the rescaled variables ω/ω c and k 2 /qL, To illustrate that this scaling is satisfied by the proposed MC algorithm, we have plotted in Fig. 11 MC simulations of ω dI dω dκ 2 as a function of κ 2 for different ranges of gluon energy, ω ∈ [ω min ; ω max ], and for different values of µ 2 (i.e. different values of q eff ), keeping L, λ el and λ inel fixed. It is then a direct consequence of (6.2) that varying ω min , ω max and q eff by the same factor will leave the distribution ω dI dω dκ 2 unchanged. Fig. 11 illustrates that this generic scaling property of the BDMPS-Z formalism is satisfied by the MC algorithm.
We finally discuss the ω-dependence of the double-differential gluon energy distribution for fixed values of κ 2 . The multiple soft scattering approximation and the N = 1 opacity approximation of (2.2) are known to result in an ω-dependence of ω dI dω dκ 2 that is flatter for increasing κ 2 [47]. The same feature is clearly seen in Fig. 12.
We finally show in Fig. 13 that the universal scaling property (6.2) is also clearly supported by the analysis of the ω-dependence of ω dI dω dκ 2 . In summary, we conclude that the MC algorithm proposed in section 4 reproduces all numerically relevant qualitative and quantitative features of the BDMPS-Z formalism.

Conclusions and Outlook
Multi-parton production processes exhibit destructive quantum interference effects. In general, their probabilistic implementation involves approximations. For multiple parton branching in the vacuum, the dominant destructive interference effect can be taken into account probabilistically by an angular ordering prescription. This probabilistic reformulation of the analytical expression is an approximation that is known to have the same parametric accuracy in log Q 2 and log 1/x as the leading order perturbative calculation. Within QCD matter, parton splitting close to the eikonal limit (2.1) is dominated by medium-induced destructive quantum interference effects that are calculated in the BDMPS-Z formalism. In this paper, we have demonstrated that the dominant medium-induced interference effect for kand ω-differential parton distributions can be taken into account probabilistically by a re-weighting of gluon emission probabilities based on gluon formation times. This probabilistic formulation is an approximation of the analytical BDMPS-Z result (2.2). We have established in a detailed numerical study to what extent it is a good approximation.
The proposed probabilistic implementation of the BDMPS-Z formalism is based on approxi-  mating by theta-functions those oscillatory functions that interpolate in the analytical BDMPS-Z formalism between the coherent and incoherent limiting cases. The proposed MC implementation reproduces by construction the known coherent and incoherent limiting cases and it interpolates, by construction, between these limits on the correct momentum scales. This ensures that the MC simulations agree in normalization and parametric dependencies with the analytically known results. In small regions of phase space and for very small in-medium path lengths, however, destructive interference effects are known to show up in the medium-induced part of the k-integrated gluon energy distribution eq. (2.2) as oscillatory behavior. The approximations in the probabilistic reformulation will not account for detailed interference effects such as oscillations in ω dI dω , but these are known to be numerically small and they depend also in analytical calculations on the approximations employed to evaluate eq. (2.2). While we have not advanced a parametric argument for the accuracy of the proposed MC algorithm, we conclude from the detailed numerical study in section 5 that the algorithm allows to implement probabilistically and in a quantitatively controlled manner all numerically relevant features of the BDMPS-Z formalism. This includes the correct normalization of the average parton energy loss and the norm and shape of the ω-differential distribution, as well as the parametric dependencies on in-medium path length, transport coefficient and coupling constant. An analogous remark applies to the k-differential distribution, as established in section 6.
The phenomenological modeling of jet quenching based on the BDMPS-Z formalism faces several longstanding problems. First, phenomenological models must account for medium-induced gluon splitting also outside the kinematic regions E ω and ω |k|, within which the BDMPS-Z formalism has been derived. Second, energy and momentum is not conserved in the BDMPS-Z formalism but its conservation at each microscopic interaction is phenomenologically important, in particular if it comes to simulating not only leading hadrons but the energy loss (a.k.a. mediummodified fragmentation) of reconstructed jets. Third, it is desirable to understand better how the medium-induced gluon radiation depends on properties of the scattering centers in the medium, and this requires the ability to vary the nature of the scattering centers in model calculations. Fourth, it is indispensable for a phenomenological model of medium-modified jet fragmentation that all components of the parton shower are treated on the same footing, and that means that all components can be subject to both elastic and inelastic interactions. In the BDMPS-Z formalism, however, radiated gluons scatter only elastically, and the projectile quark scatters only inelastically. Therefore the distribution of subleading partons obtained from the BDMPS-Z formalism must not be regarded as a suitable proxy for a medium-modified parton shower. Fifth, since the BDMPS-Z formalism has been derived close to the eikonal approximation, it is recoilness. This has resulted in a debate that distinguishes in an ad hoc way between collisional and radiative parton energy loss, rather than pushing for a physical formulation of the problem in which radiative contributions are necessarily accompanied by recoil (and therefore by effects that one typically associates with collisional energy loss).
The possibilities for improving on these major deficiencies of the BDMPS-Z formalism with refined analytical techniques appear to be limited. The proposed MC implementation of the BDMPS-Z formalism opens significant novel opportunities to this end. In particular, the proposed algorithm can be supplemented naturally with i) exact kinematics outside the region E ω |k|, ii) exact energy-momentum conservation at each interaction with the medium, iii) a large variety of models for the interaction with between projectile and medium, iv) a democratic treatment of all components of the parton shower and v) a kinematically correct, dynamical inclusion of recoil effects. In close analogy to the more mature situation in elementary particle physics, we expect that MC techniques will become in the next years also in heavy ion physics the preferred choice for the description of high-p T multi-particle final states, and we recognize their advantages in interfacing dynamical simulations of parton evolution with hadronization models. In the present paper, we have demonstrated only that the proposed MC algorithm implements all numerically relevant features of the BDMPS-Z formalism probabilistically. In our view, the main importance of this result lies in the fact that it establishes a starting point for going beyond the BDMPS-Z formalism in a framework that remains rooted in the analytically identified medium-induced interference effects. We plan to explore this approach in subsequent work.

A. Formation time of vacuum radiation from the BDMPS-Z formalism
Here, we demonstrate that a simple extension of the opacity expansion of section 2 allows one to identify within the BDMPS-Z formalism a formation time for vacuum radiation. Although this quantity does not enter the MC algorithm proposed in the present paper, we find this observation sufficiently interesting to discuss it in the present appendix.
The main idea of the following is to gain further insight into the different roles of vacuum and medium-induced radiation by introducing a length scaleL that separates the production of the partonic projectile at ξ 0 = 0 from its in-medium propagation after timeL. To this end, we study the BDMPS-Z formalism for a uniform distribution of scattering centers in a spatial region that is separated by a lengthL from ξ 0 , n(ξ) = n 0 , forL < ξ <L + L , 0 , for ξ <L or ξ >L + L . (A.1) From equation (2.2), we find then to first order in opacity the medium-induced gluon energy distribution ω dI(N = 1) dω dk dq = α s C R π 2 1 (2π) 2 |A(q)| 2 − V totδ (q) 1 (k + q) 2 + q 2 k 2 (k + q) 2 , × (n 0 L) L Q 1 − sin (L +L) Q 1 + sin L Q 1 L Q 1 . (A.2) Here, k denotes the transverse momentum of the gluon in the final state, and (k + q) can be regarded as the transverse momentum of an incoming gluonic component of the partonic projectile. The value Q 1 = (k + q) 2 /2ω denotes then the transverse energy of this initial gluonic projectile component, prior to exchanging a transverse momentum q with the medium. In the following, we investigate under which conditions this initial gluonic component can be freed (i.e. radiated) by a medium positioned betweenL andL + L. We note first that the vacuum radiation term H (k + q) in the first line of (A.2) displays the standard collinear singularity of the vacuum radiation. Also, the medium-induced radiation term R (k, q) shows singularities for vanishing incoming gluon momentum (k + q) and for vanishing outgoing gluon momentum k, as one expects for the radiation from an isolated single scattering center. We now discuss how the destructive interference term in the second line of (A.2) regulates the incoming singularity at a scale that depends on the position and thickness of the target. We consider first a medium of a fixed number of active scattering centers, that means, a medium of fixed opacity (n 0 L = fixed). For gluons of initial transverse energy Q 1 , we can then always find a sufficiently large in-medium path length L 1/Q 1 , so that these gluons can be freed with negligible destructive interference effects, L Q 1 − sin (L +L) Q 1 + sin L Q 1 L Q 1 What happens in the opposite limit, when the longitudinal extension of the medium L is small compared to the inverse transverse energy of the incoming gluon, L 1/Q 1 ? Expanding the phase factor for (L Q 1 ) 1, we find L Q 1 − sin (L +L) Q 1 + sin L Q 1 L Q 1 = 1 − cos(LQ 1 ) + 1 2 sinLQ 1 (L Q 1 ) The limit n 0 L = fixed, L → 0 corresponds to localizing medium effects exactly at a distanceL after the starting point ξ 0 of the parton evolution. In this limit, the phase factor (A.4) is 1 − cos(LQ 1 ) , and it cancels the 1/ (k + q) 2 divergencies in (A.2) only ifLQ 1 1. Therefore, gluons with initial transverse energy Q 1 can only be produced in interactions with the medium, if the medium is placed at a distanceL We note that the limit n 0 L = fixed, L → 0 can be viewed as a gedankenexperiment, according to which one produces a parton at time ξ 0 and allows for its vacuum evolution up to a timeL before testing the content of the evolved vacuum wave function by an interaction with the medium at timē L. The inequality (A.5) suggests a probabilistic picture according to which -irrespective of the nature of the medium and the strength of its interaction -one can interact with gluons of transverse energy Q 1 in the incoming vacuum wave function of the projectile only at times later than 1/Q 1 . In this sense, the inverse transverse energy 1/Q 1 of the gluonic components prior to interaction with the medium has a natural interpretation as the formation time τ (vac) f of gluons in the vacuum. Heuristic proposals for the life time of a parent parton in the vacuum are often based on its virtuality Q. In its own rest frame, a state of virtuality Q is expected to have a lifetime ∼ 1/Q. In a Lorentz frame in which this virtual partonic state has energy E, its life time ∼ 1/Q is Lorentz dilated by a boost factor E/Q, We consider now the standard perturbative situation that the virtual parent parton splits into two partons with much lower virtuality and with momentum fractions z and (1 − z) respectively. The relative transverse momentum k pair between the two daughter partons satisfies then