Parton energy loss in a hard-soft factorized approach

An energetic parton travelling through a quark-gluon plasma loses energy via occasional hard scatterings and frequent softer interactions. Whether or not these interactions admit a perturbative description, the effect of the soft interactions can be factorized and encoded in a small number of transport coefficients. In this work, we present a hard-soft factorized parton energy loss model which combines a stochastic description of soft interactions and rate-based modelling of hard scatterings. We introduce a scale to estimate the regime of validity of the stochastic description, allowing for a better understanding of the model's applicability at small and large coupling. We study the energy and fermion-number cascade of energetic partons as an application of the model.


I. INTRODUCTION
The production of energetic hadrons and jets in heavy ion collisions is markedly different from the production of energetic electroweak bosons. The latter clearly exhibit "binary scaling": weak bosons and high-energy photons are produced as if nucleons from each nucleus were independently undergoing inelastic binary collisions [1][2][3][4][5][6][7] (see also Refs. [8,9] and references therein). Hadron and jet measurements, on the other hand, display evident deviations from binary scaling. These deviations are understood to be a consequence of the formation of a quark-gluon plasma in relativistic nuclear collisions: energetic parton production does follow "binary scaling"; it is their subsequent interactions with the plasma that lead to parton energy loss, and consequently to an apparent deviation from binary scaling for hadronic observables.
This characteristic phenomena of jet and hadron "energy loss" in heavy ion collisions has been observed at both the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) [10][11][12][13][14][15]. Energetic partons are produced at the earliest stage of heavy ion collisions, and they propagate through all the different phases of the collisions. As a consequence of their interactions with the quark-gluon plasma, the momentum distribution of these energetic partons changes distinctly compared to the baseline observed in proton-proton collisions. 1 This makes them important probes of the deconfined nuclear plasma produced in heavy ion collisions.
A number of different formalisms have been used to model the interaction of energetic light partons 2 with the constituents of the plasma [17][18][19][20][21][22][23][24][25][26][27] (see also Refs. [10,[28][29][30] and references therein). Fundamentally, most parton energy loss formalisms have a well-understood common core, yet applications to heavy ion collisions tend to require approximations and practical considerations that lead to non-negligible differences between parton energy loss models [28][29][30]. One difference between the models is the treatment of the underlying plasma, which is often assumed to be made of a large number of quarks and gluons with energies of 1 GeV in near local thermal equilibrium. Whether these quarks and gluons are treated as dynamical entities or as static scattering centers is one of many differences in the energy loss formalisms [28][29][30]. The above assumption is important, given that the quark-gluon plasma produced in heavy ion collisions is understood to be strongly coupled [31], and a quasiparticle description may not be justified. 3 A different phrasing of the above challenge is that the energy loss of even very energetic partons can be affected by non-perturbative effects from the stronglycoupled plasma. Hard interactions -those with large momentum transfer between the energetic parton and the plasma -are expected to have smaller non-perturbative effects, or even be accessible perturbatively, as a consequence of the running of the QCD coupling. On the other hand, "soft" parton-plasma interactions with small momentum transfer are expected to suffer the largest nonperturbative effects. 4 We note that "hard interactions" and "soft interactions" have various meanings in the literature, but for the purpose of this work, the temperature of the plasma can be considered as the scale separating hard (larger than T ) and soft (smaller than T ) interactions.
A stochastic treatment of these soft interactions of energetic partons provides an alternative approach to account for non-perturbative effects -an approach that is agnostic to the strongly-or weakly-coupled nature of the underlying deconfined plasma. The dynamical details of the large number of soft interactions are encoded in a small number of transport coefficients. The latter can be parametrized and constrained from measurements. They can also be studied using lattice techniques (see for example Ref. [34,Section 4] and Refs. [35,36] and references therein). From a practical point of view, a stochastic description of a large number of soft interactions can also be more efficient numerically than a rate-based approach.
A systematic hard-soft factorization of parton energy loss was proposed recently to describe parton propagation in a weakly-coupled QGP [34,37]. In this factorization, soft interactions are described as a stochastic process with drag and diffusion transport coefficients calculated perturbatively; hard interactions are solved with rates that are also calculated perturbatively. In the weakly-coupled regime, parton energy loss in this hardsoft factorizated scheme was shown to be equivalent to a fully rate-based treatment of parton energy loss [37]. Importantly parton energy loss in this hard-soft factorization can also be extended to next-to-leading order [37], a feature beyond the scope of this work which we shall explore in the future.
As discussed above, the drag and diffusion contribution to parton energy loss can be factorized systematically, and calculated non-perturbatively e.g. based on Electrostatic Quantum Chromodynamics (EQCD) [35], or fitted to data. These extractions will then depend on the separation scale µ, which appears in the approach. At higher order, the drag and diffusion coefficients will evolve with the scale µ ∼ πT , incorporating in a consistent way the running of the coupling. While this is beyond the scope of this work, we hope that this manuscript can provide a first step in that direction. Throughout the paper we will already study the dependence of various observables on the separation scale µ, and, encouragingly, find that this dependence is moderate in most cases.
The above work is based on the "effective kinetic theory" approach [25] derived for a weakly-coupled quarkgluon plasma. In a weakly-coupled plasma, quark and gluon excitations are described as quasi-particles with effective properties related to the local density of the plasma. In this effective kinetic approach, the dynamics of quasi-particles are described by Boltzmann transport equations. Leading order [O(α s )] realizations of this effective kinetic approach -extrapolated to large values of strong coupling constant α s -have been used widely to study parton energy loss (see e.g. Refs. [26,[38][39][40][41]).
In this work, we present the first numerical implementation of the hard-soft factorized parton energy loss model [37] discussed above. For our implementation we utilize the publicly available JETSCAPE framework [42], as it allows us a straightforward integration of our parton energy loss model with the other ingredients necessary for a full simulation of jet production in heavy ion collisions. We first test and validate this factorization of parton energy loss in the weak coupling regime for a static medium.
We introduce a dimensionless scale to quantify the kinematic range for which soft interactions can be described accurately with a stochastic approach. We use this scale to discuss a hard-soft factorization model for a strongly-coupled quark-gluon plasma, relevant for phenomenological applications in heavy ion collisions.
Finally, we present an application of our new factorized model of parton energy loss by calculating the energy and fermion-number cascade of an energetic parton propagating in a static medium, finding good agreement with known analytical approximations. The evolution of an energetic parton in a thermal medium of temperature T can be described by a Boltzmann transport equation [43]: where P = (p, p) is the four-momentum of the energetic parton and C is the collision kernel of the parton with the medium. The index a represents partons with a certain color and helicity state. We use the same notation for the parton momentum distributions as in Ref. [37]: the distribution of rare energetic partons of type a is δf a (p, x, t), to distinguish it from the quasi-thermal distribution of soft particles n a (p, T (x, t), u(x, t)), where u is the flow velocity. In this notation, the total phase space distribution of quasiparticle a is f a (p, x, T ) = n a (p, T (x, t), u(x, t)) + δf a (p, x, t). We assume p T and g ≡ αs 4π 1. Because interactions between energetic partons themselves are rare and can be neglected, the Boltzmann equation is effectively linear in δf a .
At leading order, the interactions between quasiparticles can be divided as 2 ↔ 2 elastic interactions and 1 ↔ 2 inelastic interactions. Elastic 2 ↔ 2 processes refer to elementary scatterings involving two incoming particles and two outgoing particles without any radiation. Multiple soft 2 ↔ 2 scatterings between the energetic parton and the plasma can induce a collinear radiation. In the effective kinetic approach, these multiple soft scatterings are resummed consistently, to account for interference between subsequent collisions which lead to the Landau-Pomeranchuk-Migdal (LPM) effect. This resummed collinear radiation is known as the effective 1 ↔ 2 process. The collision kernel of both 2 ↔ 2 and 1 ↔ 2 processes can be written as Importantly, in our approach, we only follow the evolution of energetic partons with an energy above a cutoff p cut = 2 GeV. Our assumption is that we should focus our efforts on high-p T observables which are dominated by partons above this cutoff. After neglecting terms suppressed by exp(−p/T ), the collision kernels C 1↔2 and C 2↔2 read [37]: where the notation for the Lorentz-invariant integration is and ν a is the degeneracy of particle a.
For C 1↔2 , a is the incoming hard parton with the momentum p, and b, c are outgoing particles with the momentum p , k . γ a bc is the splitting kernel of a → bc, which can be calculated with the AMY integral equations [25,26]. For C 2↔2 , particle a is the incoming energetic parton with momentum p, particle b is the plasma particle with the momentum k interacting with a, and particle c, d are the outgoing particles with momentum p , k . M ab cd is the matrix element of the elementary process ab → cd [25].

B. Reformulating parton energy loss with hard-soft factorization
In the hard-soft factorized parton energy loss model introduced in Ref. [37], 1 ↔ 2 and 2 ↔ 2 processes are further divided into soft interactions and hard interactions. The collision kernel is rewritten as In this hard-soft factorized model, soft interactions described by C 1↔2 soft and C 2↔2 soft are treated stochastically with the Langevin equation.
Hard inelastic interactions, C 1↔2 hard , are treated with an emission rate as calculated from the AMY integral equations [25]. 5 We refer to them as large-ω interactions.
The hard 2 ↔ 2 part C 2↔2 hard is further divided as (i) large-angle interactions, and (ii) splitting approximation processes, based on the energy transfer ω: The physical meaning of this separation is the following. Elastic collisions occur between an energetic parton (p T ) and a lower energy quasi-thermal quark or gluon (k ∼ T ). On rare occasions, the momentum transfer in these elastic collisions is sufficient to make the low-energy Example of inelastic interaction, in which multiple soft scatterings induce the radiation of a soft gluon with energy ω. We denote the radiations with ω > µω as large-ω inelastic interactions. quark or gluon become an energetic parton with k T ; such partons are referred to in the literature as "recoil partons". The process through which a recoil parton is produced is akin to a splitting process: a single energetic particle splits in two energetic ones. The kinematic of this process simplifies and it benefits from being treated separately.
The factorization of the phase space for this reformulation is summarized in Fig. 1. An ensemble of cutoffs is used to divide the different regions of phase space. We discuss the details of the different treatments and these cutoffs in the following subsections. The diagram of a typical 1 ↔ 2 inelastic interaction is shown as Fig. 2. We assume the energy of the radiated particle is ω. We define a hard-soft cutoff µ ω based on the radiated energy ω, to divide C 1↔2 hard and C 1↔2 soft . In the weakly-coupled regime, the cutoff µ ω is limited to µ ω T , where T is the temperature of the thermal medium.
Collinear radiations with energy ω > µ ω are included into the hard part, C 1↔2 hard ; they are treated as usual with emission rates calculated from AMY's integral equations as in Eq. (3).

Treatment of hard interactions: elastic case (2 ↔ 2)
The diagram of a typical 2 ↔ 2 elastic interaction is shown in Fig. 3. We define the momentum transfer between the two incoming particles as Q = (ω, q); the four-momenta of the incoming and outgoing particles are P = (p, p), K = (k, p), P = P − Q, and K = K + Q. Usingq ⊥ ≡ q 2 − ω 2 and ω, we divide the phase space of elastic interactions as • Large-angle scattering C 2↔2 large-angle :q ⊥ > µq ⊥ and ω < Λ; • "Splitting-like" process C 2↔2 split : ω > Λ with a. Large-angle scattering (C 2↔2 large-angle ) Hard scatterings withq ⊥ > µq ⊥ and ω < Λ are denoted as large-angle interactions, because the scattering angle is generally large in this region. The cutoff µq ⊥ is typically assumed to be gT µq ⊥ T [37], although we will see in Section III that this condition can be relaxed. We assume p ω in this region, and simplify the matrix elements accordingly.
We use vacuum matrix elements for C 2↔2 hard , because the screening effects are only significant for soft interactions (C 2↔2 soft ) in the weakly-coupled regime [37]. Since we are only interested in the evolution of energetic partons, we keep terms to the first order in T /p in the matrix elements.
The treatment of the regionq ⊥ > µq ⊥ and p−Λ < ω < p -which is handled differently for technical reasonsis discussed in Appendix B.
b. Splitting approximation (C 2↔2 split ) When both of the outgoing particles of a 2 ↔ 2 interaction are hard (p , k > p cut ), the interaction can be effectively considered as a splitting process. The splitting leads to a hard recoil parton which should be included in the calculation.
We use a cutoff Λ on ω to distinguish two hard outgoing particles from only one hard outgoing particles. In principle, this cutoff Λ should be 3T Λ p. In the numerical implementation, unless specified otherwise, we choose Λ = min( √ 3pT , p cut ) to divide C 2↔2 split and C 2↔2 large−angle . Recall that we use p cut = 2 GeV in this work. As shown in Fig. 3, splitting approximation process is the 2 ↔ 2 interactions with Λ < ω < p − Λ.
At the interface between the phase space of largeangle scattering (C 2↔2 large-angle ) and splitting-like processes (C 2↔2 split ), the two collision kernels should be consistent. We verified this in Fig. 4: the differential rates of C 2↔2 split and C 2↔2 large-angle are compatible in √ 3p cut T < ω < p cut . As long as we choose the cutoff Λ in this range, this division of C 2↔2 hard should be consistent. A detailed discussion of the splitting approximation The differential rate of splitting approximation interactions and large-angle interactions for gg ↔ gg process when αs = 0.3. The shaded area is the region of √ 3pcutT < ω < pcut. dΓvac/dω is the differential rate of vacuum matrix elements for 2 ↔ 2 interactions. The results are for p0 = 100 GeV. Note that in the numerical implementation, we double-count the large-angle interaction rate because we only sample in half of the phase space. Here, to compare with splitting approximation rate, we decrease the large-angle interaction rate in the numerical implementation by a factor of 1 2 to cancel out the double-count.
process is in Appendix D. The p T and ω T kinematic cuts lead to significant simplifications for the matrix elements entering into C 2↔2 split .

Treatment of soft interactions
In the hard-soft factorized approach, the large number of soft interactions are described stochastically with drag and diffusion coefficients. When the momentum transfer is small, the Boltzmann equation [Eq. (1)] can be approximated as a Fokker-Planck equation. The collision kernel of the Fokker-Planck equation is written as: where η D,soft is the drag coefficient of the soft interactions,q L,soft andq soft are the longitudinal and transverse momentum diffusion coefficients of the soft interactions.
In the diffusion process, the number and the identity of the particles are preserved. Since the soft radiations of the 1 ↔ 2 process are absorbed by the plasma, the number of particles is also preserved. We include both the soft 1 ↔ 2 and 2 ↔ 2 collisions in the diffusion process. The diffusion process can be solved using a Langevin equation [44] in the numerical implementation.
For soft 1 ↔ 2 process, we can obtain the perturbative longitudinal diffusion coefficients by expanding C 1↔2 and only keeping the soft radiation terms. At leading order in α s ,q where C R is the Casimir factor. For gluons, C R = C A , while for quarks, C R = C F . 6 The derivation of this value can be found in Appendix A. We assume that the radiation angle is zero for collinear radiations. Consequently the transverse diffusion coefficient of 1 ↔ 2 interactions is approximated as zero.
For soft 2 ↔ 2 processes, the diffusion coefficients can be calculated perturbatively; a modern derivation can be found in Ref. [37]. The transverse momentum diffusion coefficient due to soft scatterings iŝ where is the square of the leading order Debye mass, N c = 3 is the number of colors and N f is the number of flavors involved in the interactions.
Since detailed balance is preserved in the Fokker-Planck equation, as verified in Appendix E, the drag coefficient η D can be calculated from diffusion coefficients according to Einstein relation for both soft 1 ↔ 2 and 2 ↔ 2 processes: Equations (11)(12)(13)(14) assume that the coupling α s is small. We discuss the range of validity of the perturbative coefficients in Section III A 1. Our long-term goal is to treat q soft andq L,soft as non-perturbative parameters, incorporating much more physics than leading order scattering. These parameters could then either be constrained with lattice inputs [35] or fitted to experimental data, e.g. with the Bayesian approach [46,47]. In either case, the results will depend on the separation scale µ, and this dependence would then have to match with the hard sector at LO (order g 2 ), NLO (order g 3 ), and NNLO (order g 4 , the first order the coupling runs). Ideally the hard sector, and thus the evolution with µ can be treated perturbatively. As a first step we will study the sensitivity to the scale separation µ in this manuscript.
Besides the identity preserving diffusion process, the identity of the particle can be converted through soft fermion exchange with the medium. This exchange must be screened with the non-perturbative HTL resummation scheme. In the hard-soft factorized approach adopted here, we separate the 2 ↔ 2 processes with fermion exchange into hard collisions withq ⊥ > µq ⊥ , and soft collisions withq ⊥ < µq ⊥ (see Fig. 1). The hard exchange collisions are treated with vacuum matrix elements, while the soft exchange collisions are incorporated into a conversion rate Γ conv q→g (p) for q → g: Here m 2 ∞ is the fermion asymptotic mass, m 2 ∞ = g 2 C F T 2 /4 [37,45]. In each time step there is a probability ∆t Γ conv for a quark to become a gluon, with the same momentum, and vice versa. Further details about the conversion rate C 2↔2 conv are given in Appendix C. In the future, the non-perturbative conversion coefficient Γ conv q→g can be taken from a next-to-leading order analysis [37], or can be determined from data in a Bayesian approach.

Summary
In summary, the collision kernel of hard-soft factorized model is reformulated as The cutoff dependence of the stochastic description is cancelled in Eq. (16) by the cutoff dependence of the hard interactions. That is, each individual process in the hard-soft factorized model is dependent on the cutoff, but this dependence cancels out when all the processes are summed. We show this explicitly in Section III.

C. Running of the strong coupling αs
All discussions up to this point assumed that the strong coupling constant α s is fixed at a given small value. It is clear, however, that the strong coupling constant will be different for soft and hard interactions; this is in fact a key assumption of the present model: hard interactions are more perturbative than soft ones, because the coupling constant scales inversely with the momentum exchange between the energetic parton and the plasma (see Ref. [29, Section V] for a discussion, for example). The running is slow (logarithmic in the momentum exchange), however, more studies will be necessary to understand the exact magnitude of loop corrections or nonperturbative effects on soft and hard collisions.
As a first step in introducing our model of parton energy loss, we keep the strong coupling constant α s fixed throughout the manuscript.

III. HARD-SOFT FACTORIZATION OF PARTON ENERGY LOSS IN THE WEAKLY-COUPLED REGIME: NUMERICAL STUDY
In the first part of this section, we compare the analytical equations for the soft-interaction parton transport coefficients [Eqs. (11)(12)(13)] with their numerical values evaluated from the matrix elements, and summarize the range of cutoff and coupling where they are consistent. We also compare (i) soft interactions modelled with matrix elements with (ii) soft interactions modelled with the Langevin equation. We perform this test in the weak coupling limit. We use this discussion to review the range of validity of the Fokker-Planck equation and its stochastic Langevin implementation.
In the second part of this section, we compute the energy loss of an energetic parton in a brick and discuss the dependence of the results on the soft-hard cutoffs introduced in Section II B.

A. Treatment of soft interactions
Soft interactions can be described either stochastically with transport coefficients, or microscopically with matrix elements. In what follows, we compare these two descriptions, with particular emphasis on the effect of the soft-hard cutoffs and of the coupling constant.
The tests performed in the present subsection are not expected to be related to exact composition of the plasma (number of quark flavors). Thus, for simplicity, the calculations are performed in the pure glue limit (N f = 0).

Analytical and numerical calculation of soft transport coefficients
In a weakly-coupled quark-gluon plasma, the drag and diffusion coefficients for soft interactions can be calculated analytically using perturbation theory [ Eqs. (11)(12)(13) ], as discussed in Section II B 3. The same drag and diffusion coefficients can be obtained by direct numerical integration of the parton energy loss rates; these rates are calculated from matrix elements screened by plasma effects [43].
The diffusion coefficients are defined as [37] where ∆p ⊥ is the momentum change perpendicular to the direction of the energetic parton, and ∆p L is the longitudinal momentum change of the parton. The brackets represent an average over all interactions during the parton propagation. The numerical soft diffusion rates are thus calculated as [37] q 2↔2 where dΓ(p, q)/dω and d 2 Γ(p, q)/dωdq ⊥ are the rates for an energetic parton with four-momentum (p, p) to undergo a four-momentum change (ω, q) calculated using screened matrix elements. The initial parton energy p is assumed to be much larger than all other energy scales in the problem, effectively p → ∞. The cutoffs µq ⊥ , Λ and µ ω are used to limit the phase space of interactions included in the transport coefficients, in the present case to limit the interactions to soft ones only. There are two important differences between Eq. (18) and the analytical diffusion coefficients Eqs. (11)(12)(13). First, Eq. (18) is formally valid for arbitrarily large cutoffs (µq ⊥ , Λ and µ ω ), while Eqs. (11)(12)(13) assume the cutoff to be at most of order T . Second, there is the question of the smallness of the coupling. Equations (11)(12)(13) are derived assuming α s 1. Equation (18) is valid at arbitrarily coupling, although the rates dΓ(p, q)/dω and d 2 Γ(p, q)/dωdq ⊥ themselves are typically calculated perturbatively. 7 A comparison of Eq. (18) and the analytical diffusion coefficients Eqs. (11)(12)(13) is shown in Fig. 5 as a function of the different cutoffs. This comparison is made at weak coupling (α s = 0.005) and yields the expected agreement between the two approaches, as long as the cutoffs are T . In Fig. 6 we compare the analytical soft diffusion coefficients Eqs. (11)(12)(13) with the numerical soft diffusion coefficients Eq. (18) at different values of the strong coupling 7 It is highlighted in Ref.
[29] that the AMY differential equation used to evaluate the inelastic collisions rate remains similar if interactions with the plasma are non-perturbative. One difference is the perturbative partonic collision kernel C(q) ∝ m 2 D / q 2 (q 2 + m 2 D ) that must be modified. Non-perturbative contributions to the thermal masses are another difference.  18)) and analytical (Eqs. (11)(12)(13)) momentum transport coefficients: The numerical results are computed with exact 1 ↔ 2 or 2 ↔ 2 kinematics up to a cutoff µ/T . The analytical coefficients make kinematic approximations appropriate for µ/T 1. The results are shown for different values of the hard-soft cutoffs at αs = 0.005. We calculate these results using p0 = 100 GeV and T = 300 MeV in a pure glue medium (N f = 0). The cutoff µ in the figure denotes µq ⊥ for q 2↔2 soft andq 2↔2 L,soft , and µω forq 1↔2 L,soft . In the elastic case, the additional cutoff on ω is set to Λ = min(pcut, √ 3p0T ).  6. Comparison of the numerical and analyticalq 1↔2 L,soft , q 2↔2 soft andq 2↔2 L,soft with different coupling constants αs (see Fig. 5 for description). The solid curves denote analytical results, and the circles denote numerical results. For the kinematic cutoffs, we use µq ⊥ = µω = T and Λ = min( √ 3p0T , pcut). The numerical values of the transport coefficients were calculated assuming a T = 300 MeV pure glue medium (N f = 0) and an energetic parton with p0 = 100 GeV. constant α s . We find that the analytical soft diffusion coefficients agree well with the numerical calculations even at large coupling, except for a small tension inq 2↔2 L,soft at large α s . Tension between different calculations of the soft transport coefficients are in fact not unexpected: perturbative calculations can be equivalent at order g n yet be different at order g n+1 . These differences are negligible at weak coupling, but can become significant for larger values of the coupling. This is a natural consequence of pushing the calculations beyond their regimes of validity. There is a practical consequence: two parton energy loss calculations that use the exact same approach (weakly-coupled kinetic theory) can lead to different results, when used at large coupling; neither approach is more "correct" than the other. This is important to keep in mind when comparing the present soft-hard factorized energy loss model with other implementations such as Ref. [39].

Theoretical guidance on the range of applicability of the Fokker-Planck equation
The energy loss of energetic partons through soft interactions is described by solving the Fokker-Planck equation with a stochastic Langevin approach. The applicability of the stochastic description is limited to the regime where the Fokker-Planck equation holds. This regime of applicability depends partly on properties of the interactions rates. We can summarize the regime of validity of the Fokker-Planck equation by first expanding the Boltzmann equation for soft collisions (around ω = 0): where f (p, t) is the momentum distribution of energetic partons at time t and is the k-th moment of the differential collision rate dΓ/dω. 8 By keeping only the first two terms on the right-hand side, Eq. (19) simplifies to the Fokker-Planck equation.
Assuming a single initial energetic parton of energy p 0 , the solution of the Fokker-Planck equation is The above solution simply describes the energy distribution of the energetic parton widening from scatterings withq L = ω 2 energy diffusion, and an average energy loss of ω t.
Using this solution, we can compute the ratio of the third and second terms in the expanded Boltzmann equation (Eq. (19)):  23)) around ∆p = 0, we obtain: By taking the ratio of the second and first term of this expansion, we can find the value of ∆p for which this ratio will be large: with r a constant assumed to be smaller than 1. We can use this value of ∆p as the range of momentum around the mean energy loss that can reasonably be described by the Fokker-Planck equation. Using Eq. (26) and the first term of Eq. (24), we define the scale S as When this scale S is much smaller than 1, the Fokker-Planck equation is expected to provide a good description of the Boltzmann equation in the relevant range of momentum. We emphasize that Eq. (27) was derived without any specific form for the rate dΓ/dω; in particular, the formula is the same for perturbative and nonperturbative calculations of the rate. a. Scale for inelastic rate For inelastic interactions at weak coupling, we can evaluate Eq. 27 analytically using the formula for the very soft inelastic differential rate described in Eq. (A9). In this soft inelastic limit, the scale is given by This implies that soft inelastic emissions with energy smaller than µ can be described with the Langevin equation as long as the evolution time t in the medium is sufficiently long: Assuming µ ω ∼ T results in t 1/[g 4 T ], while µ ω ∼ gT results in t 1/[gT ]. This implies that there is a very large difference between a stochastic description of soft interactions with ω T compared to soft interactions with ω gT : in the former case, one needs a plasma 1/g 3 larger. These values serve as a reminder that, while one can in principle increase the phase space of interactions described stochastically, one may need an unrealistically large plasma for this description to be valid.
b. Scale for elastic rate The dependence of the scale S [Eq. (27)] on the cutoff µq ⊥ is shown in Fig. 7, for a small and large value of the coupling constant: α s = 0.005 and 0.3. One can see that the dependence on the cutoff can be non-monotonic for small values of α s , unlike in the inelastic case. Numerical tests, as well as the analytical expression available for the second moment at small coupling [Eq. (13)], suggest that the second moment of the elastic rate is the origin of this non-monotonic dependence of the elastic scale S on µq ⊥ .

Comparison between the diffusion process and the collision rate
In this section, we verify numerically the conclusion from the previous section: we compare a stochastic and a microscopic evolution of energetic partons in a static medium. In the microscopic rate-based picture, we use kinematic cuts to forbid hard interactions of the energetic parton. Because we are comparing soft interactions, we must use screened elastic matrix elements [43] in the microscopic description. 9 The screened inelastic (1 ↔ 2) rate is obtained numerically by solving the AMY differential equation, except for very small ω values, in which case the analytical expression described in Appendix A (Eq. (A9)) is used.
We choose the hard-soft cutoffs (i.e. µ ω and µq ⊥ ) to be at the order of T in the following tests. We set the coupling to be α s = 0.005, which corresponds to g ≈ 0.25. We choose T = 300 MeV for the temperature of the fluid, and set the propagation time in the plasma to be t = (0.3/α s ) 2 = 3600 fm. 10 We perform the diffusion approach and the collision rate approach separately to calculate the single parton energy distribution of a hard 100 GeV gluon propagating in the static pure glue medium. We emphasize once again that we only include soft interactions in the test by introducing the following hard-soft cutoffs on radiation energy and momentum transfer: for C 1↔2 soft , we only include radiations with the radiation energy ω < µ ω ; while for C 2↔2 soft , we only include interactions withq ⊥ < µq ⊥ and the energy transfer ω < Λ.
According to Eq. (27), for inelastic interactions (C 1↔2 soft ) to be describable stochastically for a cutoff ∼ T , one needs t 1/[g 4 T ] ≈ 200 fm of propagation time in the conditions described above. As expected, we find in Figure 8 that for inelastic interactions, in the weakly-coupled regime, the diffusion process can reproduce the single parton energy distribution generated by the collision-rate process. The value of the scale S, shown for each cutoff µ ω , are indeed smaller than 1. As µ ω increases, small differences appear between the Langevin description and the microscopic collision approach; the scale S is correspondingly larger, though still smaller than 1.
The same results are shown for the elastic case (C 2↔2 soft ) in Fig. 9. This time, the scale S is somewhat larger, and somewhat larger differences can indeed be seen between the Langevin and collision rate descriptions. As for the elastic case, the scale S increases as the cutoff increases, where more and more collisions are described stochastically. 9 Note that this is for testing purpose only, and that this is different from the vacuum matrix elements used for C 2↔2 hard in the hard-soft factorized energy loss model. 10

B. Parton energy loss at small coupling in a static medium
Building on the validation from the previous section, we can combine our approaches for the hard and soft interactions to implement the entire hard-soft factorized parton energy loss model described in Section II B. Remember that in the following, we use vacuum matrix elements for C 2↔2 hard , since the screening effects are encoded in the drag and diffusion coefficients of soft interactions. We also extend this test to a full quark-gluon plasma, with N f = 3. We use once again α s = 0.005 (g ≈ 0.25), with a propagation time of t = (0.3/α s ) 2 = 3600 fm in a T = 300 MeV plasma.
As summarized by Eq. (16), the hard or soft processes alone are dependent on the cutoff, but their cutoff dependence cancels out when combined. We confirm that, for both the inelastic and elastic cases, the single parton energy distribution is independent on the hard-soft cut-offs at small coupling in Fig. 10, given a sufficiently long evolution time. These results are consistent with those obtained in the previous section.

IV. HARD-SOFT FACTORIZATION OF PARTON ENERGY LOSS BEYOND WEAK-COUPLING
Soft interactions between an energetic parton and a deconfined plasma are likely non-perturbative. Evaluating this non-perturbative rate from first principles is an ongoing challenge. In this section, we estimate this nonperturbative rate using a typical approach in the heavy ion literature: we use the perturbative rate and extrapolate it to large coupling.
Recall that we do not use a running coupling in this work. As such, we use the same value of α s for soft and hard interactions, with the understanding that the future introduction of a running coupling will indeed lead to smaller values of α s for hard interactions, as assumed in this work.
As discussed in Section III A 2, soft interactions can always be described stochastically, if propagation in the medium is sufficiently long. We quantified this duration as t ω 3 2 / ω 2 3 , or S 1 as defined in Eq. (27), with ω n given by Eq. (20). We emphasize once again that Eq. (27) is general, and not limited to the perturbative regime.
We can use inelastic interactions to get an estimate of the length of the medium required to describe soft interactions stochastically. When extrapolating the weaklycoupled inelastic rate to large coupling, the ω-dependence of the rate remains the same. This means that, within this approximation, the analytical expression for S -Eq. (28) -remains the same. Consequently, Eq. (29) remains the same as well, and it states that a stochastic description of inelastic interactions with ω < T requires a time t 1/[g 4 T ]. For temperatures of a few hundred MeV and a coupling g ∼ 1 − 2 encountered in heavy ion collisions, 1/[g 4 T ] < 1 fm. Under this estimate, it would be reasonable to describe stochastically soft interactions with µ T occurring in a heavy ion collision.
Note that the above conclusion is based on the estimate of the soft inelastic rate discussed above; should the non-perturbative rate differ significantly from it, it could lead to change the range of applicability of the Langevin equation. However, we do believe that the above estimates -based on extrapolations of the weakly-coupled rates to strong coupling -are encouraging.
In what follows, we use α s = 0.3 (g ≈ 2), and first compare a stochastic and a microscopic description of parton energy loss for soft interactions. We use a plasma of length 1 fm and temperature T = 300 MeV.
Note that, when the coupling is large, the analytical diffusion coefficients computed perturbatively are not necessarily consistent with numerical values obtained by direct integration of the rates (see Fig. 6 and surrounding discussion). For what follows, we use the numerical diffusion coefficients in the Langevin part of the hard-soft factorized model.

A. Comparison between diffusion process and collision rate
As in the weak-coupling case (Section III A 3), we perform this section's test in the pure glue limit (N f = 0).
We first study the inelastic interactions, and as discussed above, we expect inelastic interactions softer than ∼ T to be describable by the Langevin equation in a 1 fm brick. We show this explicitly in Fig. 11. We show calculations for three different cutoffs µ ω , and we plot the results for the scale S from Eq. (27). 11 As expected, agreement between the Langevin approach and the microscopic collision rate approach are best when S 1. In the current setting, agreement is still good for µ ω = 2T , for which S = 0.33. This is encouraging evidence that the effect of non-perturbative inelastic interactions (C 1↔2 soft ) can be treated stochastically in phenomenological applications such as heavy ion collisions.  27)). Only soft 2 ↔ 2 interactions with ω < Λ andq ⊥ < µq ⊥ are allowed. We choose Λ = min(pcut, √ 3p0T ).
The equivalent result for soft elastic interactions (C 2↔2 soft ) is shown in Fig. 12. The result is very different. On one hand, the mean energy and width of the parton distribution described with the Langevin equation is almost identical to that described with collision rates. However their shape are different, especially at smaller values of the cutoffs µq ⊥ . Agreement between the two approaches is improved when the cutoff is larger. This is also reflected in the values of the scale S, evaluated numerically with Eq. (27), which decreases with increasing µq ⊥ (see Fig. 7). This is different from what was observed (i) in the inelastic case (see Fig. 8, 11), and (ii) in the elastic case at weak coupling (see Fig. 9): both cases preferred smaller values of the cutoff. Yet this result is fully consistent with our discussion in Section III A 2 b of the scale S for the elastic rate: it is purely a consequence of the ω-dependence of the elastic rate. We verified in Fig. 13 that longer evolution times do lead to better agreement between the Langevin and the collision rate descriptions, reflected in smaller values of the scale S. Our tentative conclusion is that soft elastic collision may be more difficult to describe stochastically; it is possible that one needs a larger cutoff µq ⊥ to describe these elastic interaction stochastically, although more studies will be necessary to confirm this conclusion. Note, however, that observables which are mainly sensitive to the average energy loss and the width of the parton distribution may tolerate a wider range of soft interactions being described with the Langevin approach.
More generally, it is clear that the choice of cutoff is very important in stochastic descriptions: careful choices of cutoffs can broaden significantly the range of applicability of the factorized approach presented in this work. Importantly, the cutoff choice should be chosen based on the expected relative size of the third and second moments of the energy loss rates.

B. Parton energy loss at large coupling in a static medium
To close this section, we quantify the cutoff dependence of a 100 GeV parton propagating for 1 fm in a 300 MeV brick of plasma, with α s = 0.3. This "brick" is the same as in the previous section. The soft interactions are described with the Langevin equation, and hard interactions are included as in the full implementation of the hard-soft energy loss model (Section II B). We use N f = 3 in this test.
We plot the energy distributions with different values of the cutoff in Fig. 14. In this larger coupling regime, as expected from the results of the previous section, inelastic interactions (C 1↔2 ) are independent of the cutoff (panel (a)). For the elastic case (C 2↔2 ), the energy distributions with different values of the cutoff are slightly different in the large energy region, although the long tail of the distribution is not affected (panel (b)).
Note that we also performed a cutoff dependence test on the cutoff Λ for 2 ↔ 2 interactions. We found the energy distribution of a parton propagating in a static medium to be independent of the choice of Λ, as expected. The result and further discussion can be found in Appendix F.

V. APPLICATION: ENERGY AND FERMION-NUMBER CASCADE
In this section, we use the hard-soft factorized model to study the energy and fermion-number cascade resulting from inelastic interactions between an energetic parton and a thermal medium. This section thus focuses on C 1↔2 (Fig. 1-b) in the hard-soft factorized model; both the hard and soft inelastic interactions are included, with the soft inelastic interactions modeled by the Langevin evolution. The collision kernel C 2↔2 is switched off for this section.

A. Energy cascade of hard gluons
When a gluon propagates through a thermal QCD medium, successive medium-induced inelastic radiations result in a gluon cascade. An analytical approximation for the gluon cascade was introduced in Refs. [48,49]; it was argued that the successive medium-induced quasidemocratic emissions lead to the accumulation of gluons at zero energy and cause a power-law scaling in the small energy region. We will study this scaling in this section [50].
At leading order, the successive radiations can be assumed to be independent [51]. In the deep LPM region, where the time scale of the radiation process is much larger than the mean free path between multiple scatterings, the rate per unit time of a gluon with energy p splitting into two gluons with energy fractions z and 1−z can be approximated as 12 [51][52][53] Hereq eff is the average transverse momentum broadening of the radiated gluon, and z = ω/p with ω the energy of the radiated gluon. We have kept only the most singular parts of the splitting function at z ∼ 0. We will treatq eff as a fit parameter, and then relate it to the parameter q 2↔2 soft in Eq. (12). The energy of the initial gluon is p 0 , and we define x ≡ ω/p 0 . The evolution of the gluon spectrum D(x, τ ) =  x(dN/dx) is governed by [48,50] ∂D(x, τ ) ∂τ where and t is the evolution time of the gluon. The exact solution for Eq. (31) can be calculated via Laplace transform: As remarked in Ref. [48], this power-law gluon spectrum Eq. (33) scales as 1/ √ x in the small-x region. In order to compare with Eq. (33), we first determine the approximate value ofq eff to use in the simplified rate Eq. (30); this value also enters Eqs. (31)(32)(33). We fixq eff by comparing Eq. (30) with the full leading-order inelastic rate, as shown in Fig. 15. With parameters given in Fig. 15, we findq eff 0.04 GeV 3 at ω/p 10 −2 . We will use this value ofq eff in our analysis of the cascade below.
It should be emphasized that Eq. (30) is an approximation to the full inelastic rates corresponding to Eq. (3). Indeed, a leading-log analysis of the full rates at small z in the deep LPM regime shows that [52]: whereq 2↔2 soft is given in Eq. (12), and µ 2 ⊥ = C 0 √ 2ωq eff with C 0 ∼ 1. The cutoff µ 2 ⊥ scales with the accumulated transverse momentum of the radiated gluon over its formation time. A next-to-leading logarithmic analysis fixes the coefficient C 0 [54]: For N f = 0, p = 1 TeV, T = 300 MeV, α s = 0.1 (same as in Fig. 15), and using ω/p = 10 −2 , we can solve Eqs. (34)(35) numerically. We findq eff ≈ 0.052 GeV 3 , which as expected is close to the value we found in Fig. 15.
We next perform the gluon cascade in a pure-glue medium (N f = 0) using the hard-soft factorized model, i.e we include both soft inelastic interactions described with the Langevin equation, and rate-based hard inelastic interactions which dominate this test. In Fig. 16, we compare this numerical result calculated by the current model with the analytical spectrum in Eq. (33). We find that the numerical solution for the medium-induced cascade is reasonably well described by the approximate analytic solution. In particular, the power law behavior, dN/dx ∝ x −3/2 , is nicely captured by this solution.

B. Fermion-number cascade of gluons and quarks
The fermion-number cascade was investigated in Ref. [55]. Given the power-law scaling in the small energy region, at small x, we can write the power-law spectrum of quarks and gluons as As derived in Ref. [55], the quark-to-gluon ratio of the soft radiated partons is determined by the transformation rate between gluons and fermions. We have where K qg is the splitting function of g → qq, and K gq is the splitting function of q → gq. To test the quark-to-gluon ratio in the hard-soft factorized model, we numerically simulate the evolution of a gluon or a quark propagating through a static QGP medium (N f = 3) using the full leading order inelastic rate. We perform the calculation for both an energetic gluon and an energetic light quark with an initial energy of 10 TeV. The result is shown in Fig. 17; we find that it converges to the universal quark-to-gluon ratio when using the full QCD rates.

VI. SUMMARY AND OUTLOOK
This work introduces a new formulation of parton energy loss where soft and hard interactions with the underlying plasma are factorized and treated separately. The factorization is performed with cutoffs based on the momentum transfer of the interactions. Rare hard interactions are considered as independent successive interactions, and solved with collision rates (Sections II B 1 and II B 2); the larger momentum exchange with the medium make them more likely to be amenable to a perturbative description. On the other hand, frequent soft interactions are treated stochastically using a Langevin evolution with drag and diffusion coefficients encoding the effect of these soft interactions (Section II B 3); non-perturbative effects can thus be absorbed in these transport coefficients.
Our numerical implementation of this model (Section III) shows that this factorization works well in the weakly-coupled regime where the theory was derived [37]. In fact, by revisiting the conditions under which the Langevin equation can describe the Boltzmann equation (Section III A 2), we extended the region of phase space ("cutoffs") that can be described stochastically.
We used the dimensionless scale S (Eq. (27)) to quantify the length of a plasma necessary for soft collisions to be describable with the Langevin equation. Our numerical tests showed that this scale works very well in practice.
Because the scale S is a property of the Boltzmann equation and not a perturbative concept, we used it to extend our discussion of parton energy loss beyond the perturbative regime. We estimated that inelastic collisions resulting in parton energy loss of order T could be described stochastically in a QCD plasma of size ∼ 1 fm (Section IV). Given that inelastic interactions dominate parton energy loss for high-energy partons, this supports the applicability of the present energy loss model in heavy ion collisions.
This work paves the way to systematic phenomenological constraints on the soft transport coefficients of light partons. The key strength of our approach is that perturbative parton energy loss calculations are still being used for harder interactions -the regions of phase space where they are most likely to hold. Conversely, the interactions most sensitive to non-perturbative effectssoft interactions -are encoded in simple transport coefficients which can be constrained by comparison with measurements. A stochastic description of soft collisions can also be very efficient numerically, as a large number of soft interactions can be absorbed in the transport coefficients. These phenomenologically-constrained transport coefficients can eventually be compared with lattice results (e.g. Ref. [35]). A similar program is already being pursued for the energy loss of heavy quarks [56]; studies of light parton energy loss with a model that includes many features of soft-hard factorization, are also ongoing [47].
Future generalization of this framework includes improving the treatment of the radiation angle of collinear radiation, and the inclusion of a running coupling constant and of next-to-leading order effects; these additions will increase the type of observables that can be studied with this model. The inclusion of finite-size effects in this formalism will also be an important addition. These additions will be able to build on Ref. [47] and other works.

ACKNOWLEDGMENTS
We thank Weiyao Ke for his invaluable help in the early stage of this project, and Jacopo Ghiglieri for generously sharing notes on the splitting approximation that formed the basis of this work's discussion on the topic. We thank Sangyong Jeon, Chanwook Park, Abhijit Majumder and the other members of the JETSCAPE Collaboration for discussions regarding MARTINI and MATTER, and for their support with the JETSCAPE framework. We thank Yingru Xu for valuable discussions regarding the Langevin equation and its numerical implementation. This work was supported by the U. Appendix A: Inelastic rate at low ω At leading order, the differential rate of the 1 ↔ 2 process can be expressed using AMY's rate [26,29]: where z = ω/p and P a bc (z) are the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting kernels of the radiation a → bc, Very soft interactions (ω T ) are dominated by gluon radiations, i.e. g ↔ gg, q ↔ gq with a soft final state gluon (see footnote 6). In this case (ω T p), AMY's integral is symmetric and can be expanded in terms of the radiated energy ω [37]: where the collision kernel is .

(A4)
We define the integral in Eq. (A3) as Combining the factors, we have By rescaling all the dimensional quantities by M ∞ , the integral I can be calculated as (A7) We therefore have an analytical approximation of the soft gluon radiation rates: For the tests in Sections III and IV, we use this soft limit of differential rate when |ω| ≤ 0.2T . In Fig. (18), we compare this soft limit with AMY's full rate for g ↔ gg, and they agree well in the soft ω region.
With the soft radiation assumption ω T p, we can simplify Equation (A8) by neglecting the terms sup-pressed by ω/T and ω/p Using the above expressions, we can calculate the perturbativeq 1↔2 L, soft . We find that the longitudinal momentum broadening of soft 1 ↔ 2 iŝ Appendix B: Energy loss rate for hard 2 ↔ 2 interactions The differential energy loss rate of a hard 2 ↔ 2 interaction is calculated using the vacuum matrix elements: where ν a = 2d a is the degeneracy of the particle a,q ⊥ = q 2 − ω 2 , M ab cd is the matrix element of a vacuum 2 ↔ 2 interaction as a hard parton a interacting with a thermal particle b and transforms into particles c and d. The expression of M ab cd can be found in Table II in [25]. The collision kernel of the 2 ↔ 2 large-angle interactions is In Eq. (B3), if outgoing particles c and d are identical species, a symmetry factor of 1 2 should be included. However, this factor of 1 2 is canceled out to incorporate the interactions with p − Λ < ω < p, since symmetric 2 ↔ 2 interactions with p − Λ < ω < p are equivalent to interactions with ω < Λ. For c and d being distinct species, a factor of 1 2 is also necessary to cancel the double-count of the final states in cd . We eliminate this factor by constraining that the energy of particle c is larger than particle d. These asymmetric interactions with ω < Λ and p − Λ < ω < p are treated separately.
In Eq. (B1), the expression of Q ab cd (p, k, ω,q ⊥ , φ) is dependent on the types of particles a, b, c, and d. We summarize them using Mandelstam variables (s, t, u), Casimir factors (C A , C F ), and color degrees of freedom (d F , d A ) as follows, where C A = 3, C F = 4/3, We summarize the expression of Q ab cd (p, k, ω,q ⊥ , φ) for different interactions in Table I.
TABLE I. In this table, we use capital letter G and Q to denote hard gluons and quarks (p > pcut), and lowercase letter g and q to denote soft gluons and quarks(p < pcut). To simplify the notation, we do not specify the quark species. Q and q include the conditions of various quark species, and can also be anti-quark.
Up to order T /p, we have the following kinematics: (B4)

Appendix C: Soft conversion process
Soft conversion is a process where the identity of the hard parton is changed by its interaction with the medium. Diffusion processes only include the identity preserving soft interactions; a soft conversion process is necessary to consider identity non-preserving soft interactions.
The collision kernel of the soft conversion reads As derived in in Section 3.3 of Ref. [37], at leading or-der, the parton identity exchange rate is where m 2 ∞ ≡ g 2 C F T 2 /4 is the asymptotic mass of quarks.
Given that the rate of these identity non-preserving soft interactions is suppressed by T /p and the energy exchange ω is small, we neglect the energy loss due to these soft conversion process, and only incorporate the identity exchange.
In the numerical implementation, at each time step, we change the identity of the leading parton according to the conversion rates in Eqs. C1 and C2.

Appendix D: Splitting approximation process
As discussed in the body of the text, the collision kernel for 2 ↔ 2 scattering processes can be simplified when the energy transfer is large. 13 For simplicity, we will begin the discussion with the pure glue theory. As we will show here, and as is obvious pictorially, the 2 ↔ 2 scattering rate with large ω can be written as an effective 1 → 2 rate, which takes the form where Here we have defined and for comparison with other litterature we have defined q(δE) for the pure glue case [37] q(δE) expression we have used the thermodynamic integral, ∞ 0 dk kn B (k) = π 2 T 2 /6. The last remaining integral over z can be done and total rate for gluon absorption takes the form Γ a(g) bc where z 0 = Λ/p. The coefficients, c Λ , c p , c ln are in tab-ular form as in Table II. We note (again) that the total rate for g → gg is Γ gg /2 to account for the symmetry of the final state. We also note that the second row in this table has been summed over quark flavors.
We will now consider the case when a soft quark is absorbed from the bath, and the hard particle of type a splits a → cd. The differential rate now takes the form dΓ a(q) cd Evaluating the matrix elements (again using Table II. of [25] and Eq. (D21)), we find Again integrating over the momentum fraction we find that the total rate takes the form Γ a(q) cd = g 4 32πp where we used the integral, ∞ 0 dp p n F (p) = π 2 T 2 /12. The coefficients c Λ , c p and c ln are tabulated in Table III.   The diffusion process as described by the Fokker-Planck equation (Eq. (10)) can be stochastically realized with the Langevin model. The stochastic Langevin equations solves the evolution of the space-time coordinates and the momentum of the particle [44,57]: where x is the space coordinates of the parton, F thermal is a thermal random force satisfying the mean and the correlation function The realization of the stochastic differential equation is dependent on the discretization scheme. We choose the pre-point Ito scheme in this work [58]. In the infinite medium limit, the initial energetic partons should eventually reach the thermal equilibrium via diffusion in the thermal plasma. The equilibrium distribution of the light parton δf (p) is proportional to exp(−p/T ) in the Fokker-Planck equation (Eq. (10)), and the time derivative of the equilibrium distribution is zero. We can thus obtain the drag coefficient η D,soft as Eq. (14).
We check the thermalization of the light partons in the QGP plasma using the Langevin model (Eq. (E2)) with the drag and diffusion coefficients in Equation (11)(12)(13)(14). As shown in Figure 19, after a long evolution time, the momentum distribution of the light parton approaches the Maxwell-Jüttner distribution [59] [60]: Appendix F: Λ cutoff dependence As described in Section II, the hard elastic interactions are divided as the large-angle process and the splitting approximation process. In Fig. 20, we show the evolution of a gluon in quark-gluon plasma (N f = 3) with only 2 ↔ 2 interactions. In Fig. 20(a), with only C 2↔2 large−angle and C 2↔2 diff , the tail of the energy distribution depends significantly on the value of Λ. In Fig. 20(b), with only C 2↔2 diff and C 2↔2 split , the interactions withq ⊥ > µq ⊥ and ω < Λ is missed, which result in a missing part of the energy distribution; inevitably, the energy distribution around the initial parton energy p 0 is found to depend on Λ. In Fig. 20(c), with all the types of the 2 ↔ 2 interactions combined (C 2↔2 large−angle +C 2↔2 diff +C 2↔2 split + C 2↔2 conv ), the result is found to be independent of the cutoff Λ, as expected.