System of evolution equations for quark and gluon jet quenching with broadening

We propose a system of evolution equations that describe in-medium time-evolution of transverse-momentum-dependent quark and gluon fragmentation functions. Furthermore, we solve this system of equations using Monte Carlo methods. We then quantify the obtained solutions in terms of a few characteristic features, namely the average transverse momentum $\langle |k|\rangle$ and energy contained in a cone, which allow us to see different behaviour of quark and gluon initiated final-state radiation. In particular, the later allows us to conclude that in the gluon-initiated processes there is less energy in a cone, so that the quark jet is more collimated.


Introduction
Quantum Chromodynamics (QCD) is the well established theory of strong interactions [1].In heavy-ion collisions, nuclear matter becomes subject to extreme conditions which leads to the appearance of a quark-gluon plasma [2].This new state of matter leads to profound experimental implications, in particular for hard, or high-p T , observables such as jets [3][4][5], a phenomenon referred to as "jet quenching." In this study we plan to focus on certain aspects of jet quenching predicted in [6,7], and observed experimentally at the hadron colliders RHIC [8] and LHC [9].Jet quenching refers generally to the suppression of jets and high-p T hadrons in heavy-ion collisions due to interactions between hard partons and the quark-gluon plasma (QGP).This phenomenon is studied using various frameworks: semi-analytical [10][11][12][13][14][15][16][17], see also [18][19][20][21], and numerical methods [22][23][24] to approach parton splitting in the medium, kinetic theory [25][26][27][28][29], the AdS/CFT based approaches [30,31] and, finally, Monte Carlo methods [32][33][34][35][36][37][38][39].In terms of observables, jet quenching induces a broad range of effects that modify the internal jet substructure, contribute to out-of-cone energy loss, the decorrelation of back-to-back jets, and jet thermalization.In the recent years the particular interest was focused on the turbulent process of transferring energy from highly energetic gluon jet to soft gluons [40].This process happens without accumulation of energy by the gluons with moderate values of longitudinal momentum fraction.One of the recently actively investigated problems is the simultaneous evolution of quarks and gluons in the QGP studied in connection to turbulent behaviour introduced above.In [41] it has been demonstrated that the quark to gluon ratio of the soft fragments tends to a universal constant value that is independent of the initial conditions.In this paper we would like generalise this discussion introducing transverse momentum dependence of quarks and gluons.Besides confirming the findings of [41] our analysis will allow to have more detail information about structure of jets as well as to understand better the broadening phenomenon which is directly linked to transverse momentum dependence of fragmentation functions.In particular, in the recent study by some of us, we solved [39] the equation that takes into account momentum broadening both during branching and via elastic scattering.The resulting distribution is considerably different then the one which accounts only for broadening during elastic scattering.Furthermore, the distributions have harder spectrum than the usually used Gaussian distributions which are used to generate transverse momentum distribution factorised from distribution in longitudinal momentum.To address this problem also in the quark case, we generalise the discussion in [42,43] and obtain a system of equations linking quarks and gluons.In this approach, QGP is modelled by static centres and a jet interacts with it weakly.The jet propagating through plasma branches according to BDMPS-Z mechanism [10][11][12][13][14][15][25][26][27][28] and gets broader due to elastic scattering with plasma.Furthermore, we solve the equations in full generality, i.e. accounting for broadening during branching and due to elastic collisions.
The paper is organised as follows.In the Section 2, we derive branching kernels for quarks and gluons.The kernels allow for splitting of quarks and gluons and also account for transverse momentum dependence.In Section 3, we write a system of evolution equations for quarks and gluons and present its formal solution.In Section 4, we present distributions resulting from numerical solutions of the equations as well as results for an average transverse momentum and jet energy in a cone as a function of the cone opening angle.Then, Section 5 concludes the paper.Some further details on the evolution equations and their solutions are collected in three appendices.In Appendix A, we give explicit formulae for the splitting functions we use in our equations.In Appendix B, we provide evolution equations derived from the general ones after partial integration over some variables.For reference, in Appendix C we provide the first-order perturbative estimate of the full distributions.Finally, Appendix D contains brief descriptions of numerical methods used to solve the above equations, i.e. two Monte Carlo algorithms and a method based on the Chebyshev polynomials, as well as results of their numerical cross-checks.

Transverse-momentum-dependent splitting kernels
We shall be interested in computing the 1 → 2 in-medium splitting kernel, given by [42] K ij (Q, z, p + 0 ) = < l a t e x i t s h a 1 _ b a s e 6 4 = " h + g Z U 7 o F v g v r < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 9 q w / r 0 U T Figure 1: Illustration of the splitting function K ij (Q, z, t) for the q → g + q splitting, where where ω 0 = z(1 − z)p + 0 , and P (k) ij (z) are the (unregularised) Altarelli-Parisi splitting functions, for different QCD splitting processes. 1The three-point correlator S(3) ij in momentum space, reads where I ij refers to the path integral, Here, n(s) is the density of scattering centres in the medium, and This path integral describes the relative motion of the three internal lines of the three-point correlator during the time interval t 2 − t 1 .This can be made clear if we introduce the variables u = r i − r k and v = zr i + (1 − z)r k − r j .In this case, the effective potential takes the form where and w(q) = d 2 σel /d 2 q is the elastic scattering potential of the medium stripped of the relevant colour factor (e.g.w g (q) = N c w(q) and w q (q) = C F w(q)).In this work, we work with the thermal HTL potential, see (22) below.In Eqs. ( 5) and ( 6), the C i 's are the squared Casimir operators of the colour representation of the three correlated particles.Equation ( 6) can be proven directly by writing down the relevant 3-point functions for each individual splitting in Fig. 2.However, a more general argument relies on the properties of colour conservation [13], and can be extended to higher-order correlators as well [44].Basically, the subtraction terms in (5) correspond to the combined colour representations of the two interacting particles.For a 3-body problem, and due to colour conservation in the medium, the subtraction terms, e.g.∼ −C 0 σ(r i − r k ), necessarily have to be related to the Casimir of the particle not involved in the scattering.For a medium with constant density, i.e. n(t) = n 0 for 0 < t < L, and in the harmonic oscillator approximation, n 0 σ(r) ≈ qr 2 /4, we get an effective jet quenching parameter where Figure 2: Three of the considered colour-singlet medium correlators contributing to splitting processes in the medium.All lines represent dressed propagators resumming multiple scattering with the medium between time t 0 , corresponding to the time of splitting in the amplitude, and t 1 , corresponding to the splitting time in the complex conjugate amplitude.The two upper lines live in the amplitude, and have "mass" ω 2 = (1 − z)ω 0 and ω 1 = zω 0 , respectively from the top, while the lower line lives in the c.c. amplitude, and carries ω 0 .
Explicitly, we have It is worth keeping in mind that, in the jet quenching literature, the jet transport parameter q often refers to the gluon contribution, i.e. q ≡ qA = N c q.To summarise the results, for processes involving a gluon emission, i.e.R → g + R, where the gluon takes away the momentum fraction z, we get and in the special case of g → q + q, we get With these approximations, the solution to the path integral I ij (u 2 ; u 1 ) can be written as where Ω ij = 1+i 2 f ij (z) q/ω 0 .After performing the Fourier transforms in Eq. (2), we finally obtain which, within these approximations, is time-independent.Finally, we integrate out P and ∆t to obtain the splitting function where k 2 br = z(1 − z)p + 0 f ij (z) q is the typical transverse momentum accumulated during the time it takes to split, also called the formation time.This agrees with the expression derived in [45].
To make contact with previous works, that did not include the transverse-momentum dependence in the splitting function, we can also rewrite Eq. ( 14) as where The factor R ij (Q, k 2 br ) represents the broadening that takes place during the formation time of the splitting which is characterised by k 2 ⊥ ≈ k 2 br , and is normalised Note also that this distribution reduces to a Dirac δ-function when the transverse momentum accumulated during the formation time k 2 br tends to zero, lim Finally, in Eq. ( 15), we have also defined the analog of the stopping time for a jet with p + 0 , namely When dealing with both quark and gluon contributions to the splitting processes, this expression is stripped of the relevant colour factors.For purely gluon cascades, the correct stopping time is rather

Evolution equations and their formal Monte Carlo solutions
Using the branching kernels of the previous sections together with scattering kernels it is possible to obtain a system of equations for the evolution with time t of fragmentation functions D a for particles of type a (a = g for gluons or a = q i for quarks and antiquarks of the flavour i) or equivalently of multiplicity distributions F a , which are defined as In this section we formulate the evolution equations for the fragmentation functions and describe their Monte Carlo solutions (implemented in the MINCAS program), while the equivalent equations for multiplicity distributions have analogous Monte Carlo solutions (implemented in the TMDICE program).The evolution equations for the fragmentation functions can be written as where Q = k − zq, and the index i runs over all active quarks and antiquarks (i = 1, . . ., 2N F where N F is the number of active quark flavours in the cascade).The strong coupling constant α s in Eq. ( 20) is a function of the relative transverse momentum Q 2 ≈ k 2 br .However, we treat it here as constant: α s ≈ 0.3 (for the concrete parameter choices, see below).The elastic collision kernel C q(g) (l) is given by where is the HTL in-medium potential, where m D is the Debye mass and n the density of scattering centres in a thermal medium.At leading-order, they are given by m 2 D = (1 + n f /6)g 2 T 2 and n = m 2 D T /g 2 , where n f is the number of active flavours in the medium.The coupling to the medium g should be evaluated at the scale of the temperature ∼ 2πT .In this work, we however keep it fixed at the same value as the coupling in the medium-cascade, namely g = (4πα s ) 1/2 ≈ 2, see the concrete parameter choices below.For the HTL potential, the bare jet transport coefficient is then q = 4πα 2 s n, see e.g.[19].
In the limit of the gluon-dominated cascade, where the quark contributions can be neglected, one obtains the following evolution equation: which agrees with previous results [42].In this work, we study the interplay between quark and gluon degrees of freedom in the cascade.
In order to solve the evolution equations (20) with Markov Chain Monte Carlo (MCMC) methods, first we need to transform them into the form of integral equations of the Volterra type, similarly as in Ref. [38] for the pure gluon case.We start from introducing some useful notation that will facilitate expressing the corresponding equations in a compact and transparent form.
For transparency in the set of coupled evolution equations, let us now redefine the indices I, J to run over all parton flavours, i.e. quarks, antiquarks and gluons, so that and new Q-dependent evolution kernels K IJ (z, y, Q), defined as where 2 and K gqi ≡ K gq , K qig ≡ K qg /N F , and K qiqj ≡ δ ij K qq .The factor (1 + δ IG δ JG ) accounts for the symmetry factor that appears in the gluon-gluon splitting.The expressions for the full set of the splitting functions are given in Appendix A.
Then, let us introduce the Sudakov form-factor Ψ I (x) that resums all unresolved branchings and scatterings: where and The full branching-scattering kernel can be expressed as The analogous branching-scattering kernel GIJ of the evolution equations of the multiplicity distributions F I can be obtained with the sole replacement of With the above notation and after introducing a dimensionless evolution time τ = t/ t * , Eq. ( 20) can be cast in a simple form, namely Their formal solution in terms of the Volterra-type integral equations reads where τ 0 = t 0 / t * is the initial time for the evolution.The above integral equations can be solved numerically by iteration, where Similar solutions can be obtained for the evolution equations ( 59) and (60) given in Appendix B. In the case of Eq. ( 59) one only needs to replace in Eq. ( 30): , while for Eq.(60) in addition set w I (l) = 0.The most efficient way of numerical evaluations of the above iterative solutions is by employing the Markov Chain Monte Carlo (MCMC) methods, similar as in Ref. [38].These methods as well as the appropriate algorithms are described in Appendix D.

Numerical results
In this section we present results that we obtained as solutions of equations (20).We plot transverse 2D distributions as well as their projections.Furthermore, we use the solutions to construct characteristic features that allow us to understand better the physical effects of transverse momentum broadening and differences between quark and gluon jets.Further numerical results concerning cross-checks of different methods and programs used for solving of the above equations are presented in Appendix D. Our numerical results have been obtained for the following values of the input parameters: In the case of the initial gluon the starting distributions are: while the case of the initial quark: where the index g denotes gluons while S -the quark-singlet, i.e. the sum of active quarks and antiquarks.

Results for fragmentation functions
As a first result, we present solutions of system of equations that follows from our equations after one performs the integral over transverse momentum.The goal of this calculation is to obtain by independent calculations results presented in [41].The integrated fragmentation functions reads Since we present the results in terms of k T ≡ |k| dependence, we introduce the distribution To solve the system of the evolution equations we have used two Monte Carlo programs, MINCAS and TMDICE, as well as the numerical method based on the Chebyshev polynomials.The corresponding algorithms are explicitly presented in Appendix D. The results of Monte Carlo programs agree very well, therefore here we present the results obtained only by MINCAS.The plots comparing the Monte Carlo solutions obtained by MINCAS and TMDICE are shown in the Appendix D. Furthermore, we compare the distributions to a perturbative estimate in Appendix C. Furthermore, by visually comparing our results in Fig. 3 to the ones presented in [41], we see that we get the same features of the distributions, i.e. energy is not accumulated at the moderate values of x.Simultaneously, the distributions increase at small x values, following roughly the 1/ √ x behaviour for the gluons and the quark-singlet for the case of the initial gluon and a similar behaviour for the case of the initial quark.In the case of the initial gluon at late times there is a region, at high x, where quarks dominate.In the case of the initial quark, gluons tend to dominate if x is low for all time scales, while quarks dominate at x > 0.5.In Fig. 4, we show the k T -dependent distributions integrated over the longitudinal momenta, given in Eq. ( 40).One can clearly see that as the time progresses the distribution for both quarks and gluons become wider.Furthermore, the distributions of gluons are higher than that of quarks if gluons are in the initial state, and similarly, the distributions of quarks are higher than that of gluons if quarks are in the initial state.
The complete 2D distributions, visualising both the x and k T -dependence, are presented in Fig. 5 and 6 for the initial gluon and quark, respectively.One can can see that the late time behaviour is very similar in processes initiated by quarks or gluons.This behaviour can be linked to diffusive properties of the jet-medium interactions.Furthermore, while there are some differences in large x part of the spectra, the shape at low x is rather universal.

Multiplicity distributions
We have also obtained the results for specific times t in terms of the multiplicity distributions for jets in the medium with the TMDICE algorithm.For the numerical calculations, the same constraints as given in Eqs.(35), (36), and (37) were used, with the exception of x min , where x min = 10 −2 was chosen.The multiplicity distributions are not infrared and collinearly safe observables, and they considerably depend on the cut-off scale x min .In the solution to the evolution equations ( 20), x min corresponds to an energy scale at which the assumption that coherent medium-induced radiation and scatterings occur and dominate breaks down.Thus, we assume for this energy scale an estimate xE = 1 GeV which is still larger than the medium temperature (where we assume that jet-particles thermalize), but of the same order.The numerical results are shown in Fig. 7 for the time scales of t = 2 fm and 4 fm, for jets initiated by either a quark or a gluon.As can be seen for the distribution in x, the peak in the infrared region that is suppressed by a factor x in the fragmentation functions D(x) in Fig. 14 is much more pronounced in Fig. 7, showing that the infrared contributions to the jet-multiplicity dominate at large time scales.The distributions in k T exhibit an increasing broadening with large time scales that is more pronounced for the jets initiated by gluons than those initiated by quarks.

Characteristic features of cascades
In this section, we present our results for some useful characteristic features of the cascades.These should not be confused with observables that can be measured directly in experiment, but rather as projections or "summary statistics" containing the most important information of the cascades described in the previous section.
The first feature we discuss, is the average transverse momentum.It is defined as where k T ≡ |k|, as a function of x evaluated for different evolution time values.Overall, we see in Fig. 8 that the distributions both for the initial-state gluons and initial-state quarks are rather similar, i.e. as x gets smaller and smaller the average k T gets smaller, meaning that soft mini-jets become delocalised.The distributions for different times tend to merge as x gets small enough.One can see certain differences if one compares the distributions for the final quark vs. the final gluon as time progresses.The slopes of the distributions become different.This is   consistent with the obtained earlier results for the fragmentation functions and suggests that gluons have harder momenta than quarks and dominate at larger values of x.
In order to draw further attention to small angles, close to the direction of the parent parton, we plot the average angle θ as a function of the momentum fraction x in Fig. 9, where θ (following the light-cone kinematics) is defined as This sheds light on the internal structure of jets.We observe that gluons in the cascade typically occupy larger angles than the quarks.The average angle for gluons also grows as we go to smaller x, which is caused by broadening.This is in line with earlier observations.We also note that the average angles are very similar for gluon-initiated and quark-initiated cascades from early times.The next characteristic feature that we study is the energy contained inside a cone-angle Θ around an initial parton defined as The observable measures the amount of energy that is contained in a cone.Clearly, we see in Fig. 10 that the configuration that maximises this observable is the one where the type of the parton does not change.Furthermore, since the quark contains more in-cone energy, we conclude that it is more collimated.We can compare this to an analytical Ansatz of how the distribution should look like: where we approximate the x distribution with the analytical solution for a purely gluonic cascade [40], This Ansatz should, in principle, only be compared to the gluon distribution resulting from a fragmenting gluon.
Besides this, we compare the results to where P (k, t) is given by Eq. ( 49), while the x distribution is given by a numerical solution of Eq. (60), i.e. the evolution equations for the energy distribution.All of the results for the initial gluon feature universal slow growth if the angle is large enough, while for the processes initiated by the quark we see that the energy saturates.Furthermore, by comparing the fully analytical Ansatz D A (x, k, t) with other results, we see that both agree at large angles, while at small angles, the analytical results overestimate the amount of energy in the cone.The fully analytical Ansatz D A (x, k, t) predicts more collimated jets than the other distributions.This is because, as realised in [38], the analytical distribution given by the Gaussian function P (k, t) is narrower than the one obtained from the exact numerical solution which can be viewed as superposition of many Gaussians with different widths.The above plots allow us also to conclude that the k T dependence controls the angle at which the distribution starts to saturate.
The final feature that we discuss is the normalised amount of energy in a cone as a function of time.The results are presented in Fig. 11.One can clearly see that for processes with quarks in the initial state, quarks dominate at any time-scales in both scenarios for the angles that we consider.But what is more striking is that even when gluons are in the initial state, quarks start to dominate at late times.One can also see that the analytical solution Eq. ( 47) overshoots the Monte Carlo solution of Eq. ( 23) for the short time-scales while underestimates it for the long time-scales.

Conclusions
Our goal in this work was to study simultaneous evolution of quarks and gluons with transverse and longitudinal fragmentation functions in medium.In order to achieve that, first of all, we have introduced transverse-momentumdependent splitting kernels that take into account broadening during branching.The derivation assumes that there is a local in time factorisation between longitudinal and transverse momentum dependence, valid when the transverse momentum is much smaller than the total energy of the collision.The obtained system of equations has been solved in full generality using Monte Carlo methods as well as by the Chebyshev-polynomial method applied to the equations for the energy distribution.
] (solid lines), the Gaussian approximation in Eq. ( 50) (dotted lines) and the analytical Ansatz in Eq. ( 47) (green dashed line); the red lines correspond to the final gluon jets and the blue lines to the final quark jets.The curves for the analytical Ansatz are plotted only in the case of pure gluon jets.
Using the obtained solutions we have defined some characteristic features that allow us to demonstrate that • as evolution in time progresses gluons broaden more than quarks, see Fig. 4; • at late times, for both quarks and gluons, there is universal distribution in k T at low x; • quarks are more collimated than gluons -this we conclude from the calculation of energy contained in a cone (in the case of quarks more energy is contained in a cone); • quarks dominate at late times, see Fig. 11.
We should add here that our study breaks down at x ∼ 10 −5 where one should account for thermalisation and possible recombination of gluons, solving an appropriate Boltzmann equation that also accounts for elastic rescattering.
As an outlook, we mention here several improvements that are natural to consider for future work.A natural step is to account for medium expansion.This was already implemented for a purely radiative cascade in [46,47], which made use of appropriate splitting kernels in an expanding background.Furthermore, it would be interesting to extend the formalism presented in this paper to include hard emissions, generated by rare, hard scatterings, as developed in [18][19][20][21]48] which also is well suited for expanding media, see also [23,24,49,50] for numerical approaches to evaluating the splitting kernel.Finally, we also notice that coherence effects are important to account for finitesize medium corrections [51].These extensions are, however, beyond the scope of the current paper.

A Splitting functions
Here we provide exact equations for the splitting functions we are using.The main object is the splitting function where and where the functions f IJ (z) can be read off from Eq. ( 9).Finally, the z-dependent in-medium splitting functions are defined as where we use the following expressions for the unregularised Altarelli-Parisi splitting functions:

B Other forms of evolution equations
For completeness, we present here also the equations obtained from Eq. ( 20) by neglecting the transverse momentum accumulated during branching, k 2 br ≈ 0, and integrating over the momentum transfer Q, we get where τ = t/ t * and the sum over i runs over both quarks and antiquarks.By integrating out completely the transverse momentum, we obtain We can also rewrite the latter equations in the gluon/singlet/non-singlet basis [41], where we find with At this level, we see that it is quite natural to absorb an extra factor 2 into the gluon-gluon and gluon-quark splitting functions, respectively.Compared to the notation in [41], we identify their splitting functions with ours as: K gg (z) = 2P gg (z), K qg (z) = 2P qg (z), K gq (z) = P gq (z) and K qq (z) = P qq (z).

C First-order perturbative estimate of the distribution functions
It is very instructive to compare the result of the full medium evolution with a perturbative estimate which only accounts for one single emission or elastic scattering with the medium.In this Appendix, we summarise the results for such a first-order perturbative estimate based on formulas .
For the partonic gluon initiator, we find the first-order perturbative estimate to be, where . The term in curly brackets next to the initial condition D (0) g (x, k) is responsible for probability conservation.Now, integrating out the transverse momentum k or the momentum fraction x, respectively, we get the following distributions: which respects 1 0 dx D(x, t) = 1 and d 2 k D(k, t) = 1.The expressions for these take the form, and, finally, For readability, we have dropped terms that are proportional to the initial condition.They are however crucial in order to restore the normalisation of the distributions.Finally, in order to make contact with Eq. ( 41), we note that D(k, δt) = D(k 2 , δt), so that D(1) where k T ≡ |k|.For the quark initiator, we obtain, The integrated distributions now read, omitting the terms proportional to the initial condition, D (1)  q (x, δt) = δt α s xK qq (x, p + 0 ) , and, finally, D Numerical-solution methods

D.1 Monte Carlo algorithms
Below we briefly sketch appropriate algorithms implemented in two independent Monte Carlo programs, MINCAS and TMDICE.One essential difference between these approaches is that MINCAS provides Monte Carlo solutions for the fragmentation functions that follow Eqs.(20), while TMDICE solves the equivalent equations for the multiplicity distributions, i.e. with the replacement: First, let us rewrite Eq. ( 33) in a probabilistic form which is suited for the MCMC algorithm.To this end we define two probability distribution functions (p.d.f.s): for generation of the random variable τ i , and for generation of the random variables (z i , Q i , l i ), and the conditional probability for generating a new parton flavour J i , given the flavour J i−1 : One may easily check that they are properly normalised to 1.Then, Eq. ( 33) can be written as where the initial functions are The formula of Eq. ( 82) is a basis of the following MCMC algorithm: 1. Set the initial values of (x 0 , k 0 ) and J 0 or generate them according to the probabilities χ(x 0 , k 0 ) and p J0 , respectively.
2. Having generated the random-walk leap i − 1 (i = 1, 2, . ..), generate τ i according to the p.d.f.Ji−1 (τ i ): if 3. Otherwise, i.e. τ i ≤ τ : (a) generate the parton flavour J i according the probability p Ji , (b) generate the variables (z i , Q i , l i ) according to the p.d.f.Ξ increment i → i + 1 and go to step 2.
For the Monte Carlo solution of the evolution equations for the multiplicity distributions, as implemented in TMDICE, the algorithm is analogous, while the probability distributions change: instead from Ξ JiJi−1 (z i , Q i , l i ), the variables in step 3(b) of the algorithm are generated from the distribution where GIJ (z, y, where KIJ (z, y, Q) = K IJ (z, y, Q)/z.Please note that the functions Φ JiJi−1 , Ψ JiJi−1 , and W Ji are the same in both the Monte Carlo solution methods, TMDICE and MINCAS.
The Monte Carlo weight corresponding to the n-leap trajectory γ n is One can prove that the expectation value of such a weight corresponds to the distribution function D I (x, k, τ ): In Monte Carlo computations, based on the law of large numbers (LLN), the above expectation value is approximated by the arithmetic mean of the Monte Carlo weights for the given (x, k, τ ).In practice, one makes histograms of desired variables with the event weight given by Eq. (86).In the above MCMC algorithm, it is assumed that all the necessary random variable can be sampled from the respective probability distribution functions.In practice, however, this is not straightforward, as the corresponding evolution kernels used to construct these p.d.f.s are complicated functions, in particular they cannot be integrated analytically.Therefore, one cannot apply basic Monte Carlo techniques, such as analytical inverse transform method, to generate the respective random variables.One possible solution is to perform numerical integration and use numerical inverse transform method for random variables generation.This, however, can be slow in terms of CPU and prone to numerical instabilities, therefore it has to be done with great care.These methods have been implemented in the Monte Carlo program TMDICE for the case of the the evolution of the multiplicity distributions.
The other way to deal with this problem is to replace the exact p.d.f.s in Eq. (82) with their approximations, i.e. some simpler functions that on the one hand are as close as possible to the original functions, but on the other hand can be integrated analytically and used in the analytical inverse transform method for the random variables generation.In such a case, these approximate p.d.f.s are used to generate all the necessary random variables, and all the simplifications are compensated by appropriate Monte Carlo weights being the ratios of the exact to approximate p.d.f.s.Here, two Monte Carlo algorithms are possible: (1) the weighted-event (or variable-weight) algorithm in which each generated event is associated with the corresponding (variable) Monte Carlo weight and (2) the veto algorithm [53] which uses a dedicated acceptance-complement method to generate unweighted (i.e.weight = 1) events.These two algorithms have been implemented in the Monte Carlo event generator MINCAS.They have comparable efficiency in generation of (x, k) distributions in terms of the CPU time.It has been checked that both these algorithms give the same numerical results, which constitutes an important internal test of their implementation in MINCAS.
More details on the above Monte Carlo algorithms as well as the programs MINCAS and TMDICE will be given in separate publications.

D.2 Chebyshev method
As a test of the obtained Monte Carlo solutions, we have used semi-analytical method based on the Chebyshev polynomials to solve equations (61).The idea is to expand the solution with the Chebyshev polynomials: where • T i are the Chebyshev polynomials of the first kind, • means that the first term in the sum is divided by 2, • {y i } i∈ 1;N are the nodes (or zeros) of T N : 2 ), i ∈ 1; N , • in the following, we use the notation: The idea is then to solve the equations in the nodes of the Chebyshev polynomials in τ = t/ t * , N S (x k , τ ) ∂τ The integrals from 0 to 1 in (61) are split in integral from 0 to x k and from x k to 1 and evaluated with another Chebyshev decomposition (with adequate bijections).Actually, we use the logarithmic bijections (implying a low-x cut-off) which are more suitable here.Finally, a simple Euler method is used to solve the equation in time, with a narrow Gaussian function as a starting distribution: where ε = 10 −2 , as the Chebyshev polynomials are unable to reproduce the Dirac δ-function.The inability to approach peaked functions translates into oscillations of the solution in the low evolution times (particularly visible on the quark distributions for the quark-initiated cascade).
Unfortunately, this method has not shown suitable for the k T -dependent evolution equations (20).

D.3 Numerical cross-checks
Here, we present some results of numerical cross-check of the presented above methods for solving of the evolution equations ( 20) and (61).They have been obtained for the input parameters given in Section 4.  (20) with the results from the Chebyshev method for the evolution equations (61) at the timescales t = 0.1 and 4 fm: cascades initiated by gluons (left) and quarks (right).The bottom panels show the ratios to the corresponding results from MINCAS.
We start from the D(x, t) distributions for which we can compare all three numerical-solution methods.In Fig. 14 we present comparisons for the results from the Monte Carlo programs MINCAS and TMDICE for the evolution equations (20) with the results from the Chebyshev method for the evolution equations (61) at the timescales t = 0.1 and 4 fm.The results from MINCAS and TMDICE are obtained for the full (x, k) evolution but they are integrated by Monte Carlo method over k.One can see a very good agreement of the three methods for both the final gluon and quark distributions in both the gluon and quark initiated cascades.In the case of the Monte Carlo programs, this also shows that they reproduce correctly the D(x, t) distribution from the full evolution for the D(x, k, t) distribution given by Eq. (20).A similar agreement has been found for other values of the evolution time values.The Chebyshev method reveals some oscillatory behaviour for t = 0.1 fm, i.e. for short evolution time, but works well for higher time-scales, t > 1 fm.
In Fig. 15, we present comparisons of the D(k T , t) distributions for the evolution equations (20) obtained from two Monte Carlo programs: MINCAS and TMDICE, for the cascades initiated by gluons (left) and quarks (right) at evolution time-scales t = 0.1 and 4 fm.As one can see, a perfect agreement between the two independent Monte Carlo algorithms and programs is found.We have checked that such an agreement holds also for other evolution timescales, other splitting kernels and other medium-scattering functions.Unfortunately, for the equations including

1 k 5 w b 5 C
T 2 P I 5 Z i s = " >A A A C 7 H i c b Z L d b t M w F M f d 8 D X K V w u X 3 E R U S F x V C Z o G l x M I C Y m b V b T b t C a K H O d k N X M c y z 4 B K i t v w R 3 a 7 X g Q n o K 3 w U m T i m 4 c K d J f v / P h 4 7 + T K s E N B s G f g X f r 9 p 2 7 9 / b u D x 8 8 f P T 4 y W j 8 9 N i U l W a w Y K U o 9 W l K D Q g u Y Y E c B Z w q D b R I B Z y S Z I 2 5 p W i e P b z + y D f F 8 e t p e D D d n + 1 P D t 9 1 P 8 A e e U 5 e k F c k J G / I I f l I j s i C M P K N X J I r 8 s u T 3 g / v p 3 e 5 K f U G X c 8 z s h P e 1 V 8 F B P O x < / l a t e x i t > K ij < l a t e x i t s h a 1 _ b a s e 6 4 = " J A 1 k c n 3 R j P 1 4 U Y 4 d l D 4 y r f R M w g M = " > A A A C 2 n i c b Z J L a x s x E M f l 7 S O J + 0 j S H n t Z a g o 9 m d 0 S 2 h x D e u k x B j s J 8 S 5 G K 8 3 G I l p J S L O l R u y l t 5 J r + o n 6 K f p t K t u 7 J k 4 6 I P T n N w

Figure 3 :
Figure 3: The √ xD(x, t) distributions at the time-scales t = 0.1, 1, 4 fm: cascades initiated by gluon (left) and quark (right).The dashed lines correspond to the quark distributions while the solid lines to the gluon distributions.

Figure 4 :
Figure 4: The D(k T , t) distributions for w(l) ∝ 1/[l 2 (m 2 D + l 2 )] at the time-scales t = 0.1, 1, 4 fm: cascades initiated by gluon (left) and quark (right).The dashed lines correspond to the quark distributions while the solid lines to the gluon distributions.

Figure 11 :
Figure 11: Time evolution of the jet energy in a fixed cone with θ = 0.1 (upper panels) and θ = 1 (lower panels) for the initial gluon jet (left panels) and the initial quark jet (right panels) in the case of the k T -dependent branching and w(l) ∝ 1/[l 2 (m 2 D + l 2 )] (solid lines), the Gaussian approximation in Eq. (50) (dotted lines) and the analytical Ansatz in Eq. (47) (green dashed line); the red lines correspond to the final gluon jets and the blue lines to the final quark jets.The curves for the analytical Ansatz are plotted only in the case of pure gluon jets.

Figure 12 :
Figure 12: The √ xD(x, t) distributions at the time-scales t = 0.1, 1, 4 fm: cascades initiated by gluon.The solid lines correspond to the full medium evolution and the dotted line to the perturbative estimate.

Figure 13 :
Figure 13: The D(k T , t) distributions at the time-scales t = 0.1, 1, 4 fm: cascades initiated by gluon.The solid lines correspond to the full medium evolution and the dotted line to the perturbative estimate.

Figure 14 :
Figure 14: Comparisons of the x distributions from the Monte Carlo programs MINCAS and TMDICE for the evolution equations(20) with the results from the Chebyshev method for the evolution equations (61) at the timescales t = 0.1 and 4 fm: cascades initiated by gluons (left) and quarks (right).The bottom panels show the ratios to the corresponding results from MINCAS.

Figure 15 :
Figure 15: Comparisons of the k T distributions from MINCAS and TMDICE for the evolution equations (20) with w(l) ∝ 1/[l 2 (m 2 D + l 2 )] at the time-scales t = 0.1 and 4 fm: cascades initiated by gluons (left) and quarks (right).The bottom panels show the ratios to the corresponding results from MINCAS.