Solutions of evolution equations for medium-induced QCD cascades

In this paper we present solutions of evolution equations for inclusive distribution of gluons as produced by jet traversing quark–gluon plasma. We reformulate the original equations in such a form that virtual and unresolved-real emissions as well as unresolved collisions with medium are resummed in a Sudakov-type form factor. The resulting integral equations are then solved most efficiently with use of newly developed Markov Chain Monte Carlo algorithms implemented in a dedicated program called MINCAS. Their results for a gluon energy density are compared with an analytical solution and a differential numerical method. Some results for gluon transverse-momentum distributions are also presented. They exhibit interesting patterns not discussed so far in the literature, in particular a departure from the Gaussian behaviour – which does not happen in approximate analytical solutions.


Introduction
Quantum chromodynamics (QCD) is the well established theory of strong interactions. However, there are QCD phenomena that still require better understanding. One of such phenomena is jet quenching predicted in [1,2] and already observed in the context of the RHIC physics [3] (for an overview see [4,5] and references therein), i.e. stopping of a hadronic jet produced in an early stage of heavy ion collisions and propagating through quark-gluon plasma (QGP) which is formed in a later stage of the collisions. With the LHC being in operation, the jet quenching can be observed at much higher available energies in collisions of lead nuclei [6]. the jet-QGP interaction mechanism and a pattern of energy loss. Various approaches have been proposed which differ in assumptions about properties of plasma and jet-plasma interactions. Examples are: the kinetic theory assuming that the jet-plasma interactions can be described within a weakcoupling regime of QCD [7][8][9][10][11][12][13][14][15][16], the AdS/CFT models where one assumes the plasma to be strongly coupled [17,18] or the classical-field-theory-based approach [19] (for reviews see [7,[20][21][22][23]). Some of the mentioned formalisms are implemented in Monte Carlo event generators [24][25][26][27][28][29].
In this paper we look closer at the results obtained in [30,31] and focus on an analysis of the generation of transverse momenta via cascades of subsequently emitted jets from an energetic jet traversing QGP. In this approach the plasma is modelled by static centres and the jet interacts with it weakly. Using equations for the energy loss of the jet traversing QGP, the authors of [30][31][32] found the process to have turbulent properties, i.e. the energy is transported from large values of x to low values of x without being accumulated at intermediate values. In this paper we investigate a more exclusive equation, i.e. the equation which describes time evolution of longitudinal as well as transverse momenta distributions of gluons emitted from the energetic jet. So far this equation has not been solved numerically, and the analytical as well as numerical analyses are limited to some special cases [31,33,34] where for instance part the equation leading to broadening of transverse momentum is simplified or it is included as an input distribution. The simplified analysis suggests that the distribution of gluon transverse momenta being a solution of the equation remains Gaussian [33]. We, however, find that the exact numerical solution of the evolution equation given in Ref. [31] is not Gaussian. Furthermore, the numerical method allows to test consequences of the assumptions about properties of the medium for distribution patterns of the jets emitted from the hard jet. The paper is organised as follows. In Sect. 2, we introduce and overview the equations for the jet energy distribution and for the inclusive gluon distribution. In Sect. 3, we present reformulation of the above equations making use of Sudakovlike resummation, i.e. we resum virtual and unresolved-real emissions as well as unresolved collisions with the medium of minijets from the highly energetic jet in form of a Sudakovtype form factor. Then, we provide formal iterative solutions of these equations. In Sect. 4, we propose a Markov Chain Monte Carlo (MCMC) algorithms for numerical solutions of the above equations. In Sect. 5, we describe a numerical algorithm for solving the integro-differential equation for the jet energy distribution which is based on application of the Runge-Kutta method and discuss its limitations in obtaining high accuracy solutions. In Sect. 6, first, we present numerical results from the MCMC method for the jet energy distribution and compare them with an analytical solution as well as with results from the differential Runge-Kutta-based method. Then, we show and discuss some results for the jet transverse-momentum distributions obtained with the MCMC method. We summarise our work and present its outlook in Sect. 7. Finally, in Appendix A we provide some further details on the MCMC algorithms, in particular we describe a combination of the branching Monte Carlo method with the importance sampling.

Evolution equations
The evolution equation for gluon transverse-momentumdependent distribution D(x, k, t) in the dense medium, obtained under the assumption that the momentum transfer in the kernel is small, reads [31] where is the z-kernel function, and where t * is a stopping time, i.e. the time at which the energy of an incoming parton has been radiated in form of soft gluons, E is the energy of the incoming parton, x -its longi-tudinal momentum fraction, k = (k x , k y ) -its transversemomentum vector,q -the quenching parameter, α s -the QCD coupling constant and N c -the number of colours.
The collision kernel C(q) is given by where the function w(q), which models out-of-equilibrium momentum distributions of medium quasi-particles, takes the form [31] with q = (q x , q y ) being transverse-momentum vector and n -the density of scatterers. However, we can also consider a situation where the quark-gluon plasma equilibrates and the above transverse-momentum distribution assumes the form [35] w where m D is the Debye mass of the medium quasi-particles. In the following we shall consider both the above expressions for w(q). After integration of Eq. (1) over the transverse momentum k one obtains the evolution equation for gluon energy density [31] This integro-differential equation can be solved analytically for a simplified case of f (z) = 1 and D(x, t = 0) = δ(1−x) [31]: where τ = t/t * .

Integral equations and iterative solutions
Let us rewrite Eq. (1) by moving all terms with the minus sign from RHS to LHS and using τ = t/t * : Then, after introducing the following notation: we can write where we have introduced the upper (infra-red) cut-off for the z integral and the lower cut-off q min in the integral of w(q) over q, since the former is divergent for z → 1 and the latter is divergent for |q| → 0 in both cases of w(q) given in Eqs. (5) and (6). The whole equation does not depend on these cut-offs, so they (if sufficiently small) only play roles of (dummy) regulators of the corresponding integrals. The expression for the integral W depends on the actual form of the function w(q). For w(q) given in Eq. (5) we get while for the one of Eq. (6): The latter features much milder dependence on the lower cutoff q min than the former. The integral (x) does not have a compact analytical form, instead it can be computed numerically for a given value of . However, as will be seen later on, we do not need its actual value. The above integro-differential equation can then be transformed into an integral equation where we have introduced the following notation The factor e − (x)(τ −τ ) is the Sudakov-type form factor corresponding to resummation of virtual and unresolved-real gluon emissions as well as unresolved collisions with the medium due to the kernel-function C(q). Similarly, Eq. (7) can be transformed into the integral equation In this case the Sudakov-type form factor e − (x)(τ −τ ) corresponds only to resummation of virtual and unresolved-real gluon emissions, as there is no medium-collision term. The above integral equations can be formally solved by iteration. For Eq. (16) we obtain where with x 0 and k 0 being some initial values of x and k at the initial evolution time τ 0 , given by the distribution D(x 0 , k 0 , τ 0 ).
A similar solution can be found for Eq. (18):

Markov Chain Monte Carlo algorithms
The formal solutions given in Eqs. (19) and (21) can be used to develop Markov Chain Monte Carlo algorithms for numerical evaluation of the distribution functions D(x, τ ) and D(x, k, τ ), given some initial functions D(x 0 , τ 0 ) and D(x 0 , k 0 , τ 0 ), respectively. In Eqs. (19) and (21) there is ordering in the variable τ i : Therefore, this variable can be treated as an evolution time of a random walk in the MCMC algorithm to be used for a numerical solution of the respective integral equation. This random walk will start at some time τ 0 and finish at the time moment τ , making an arbitrary number of random leaps between τ i−1 and τ i , i = 1, 2, . . . In order to construct such an algorithm let us rewrite this equation in a probabilistic form.
First we notice that and this can be regarded the probability of a single jump beyond τ from the point τ n , i.e. a stopping rule for the random walk. Thus, the probability density function (pdf) of a random variable τ i for a single leap from τ i−1 is The random variable τ i can be generated according to the above pdf using the analytical inverse transform method.
The pdf for the variables z i and q i is given by The variables z i and q i can be generated almost independently using the following branching MC method: • with the probability p = ( and set q i = 0, • otherwise, i.e. with the probability 1 − p, set z i = 1 and generate q i according to the density function While the random variable q i can be generated according to the pdf ω(q i ) for w(q) given in Eqs. (5) and (6), generating z i according to the pdf ζ(z i ) is more difficult due to the complicated function K(z) given in Eq. (2). For this purpose one can use the rejection method or the importance samplingwe shall come back to this later on. Finally, let us define the probability density for the initial variables x 0 and k 0 : If this function is complicated, for generation of the random variables x 0 and k 0 one can use some self-adaptive Monte Carlo (MC) sampler, e.g. FOAM [36]. However, quite often it factorises into a product of probability densities: where for υ(k 0 ) one can use e.g. the Gaussian distribution: which can be easily generated, e.g. using the Box-Muller method.
Having defined all the necessary probability distribution functions, we can rewrite Eq. (19) in the following form Now we can propose the following MCMC algorithm for numerical numerical evaluation of Eq. (31): Step 1 Start a random walk from the point τ 0 . First generate the variables x 0 ∈ [0, 1] and k 0 according to the probability density η(x 0 , k 0 ), then generate τ 1 ∈ [τ 0 , +∞) according to the probability density ρ(τ 1 ). If τ 1 > τ, set x = x 0 , k = k 0 and stop the random walk, otherwise go to step 2.
Step 2 Generate the variables z 1 ∈ [0, 1] and q 1 > q min according to the probability density ξ(z 1 , q 1 ) and calculate 1 and stop the random walk, otherwise go to step 3. ...
Step n Generate the variables z n ∈ [0, 1] and q n > q min according to the probability density ξ(z n , q n ) and calculate x n = z n x n−1 , k n = q n + z n k n−1 . Then generate τ n+1 ∈ [τ n , +∞) according to the probability density ρ(τ n+1 ): if τ n+1 > τ, set x = x n , k = k n and stop the random walk, otherwise go to step n + 1. ...
Repeat the above steps N -times histogramming the variables x and k. At the end normalise the histograms with the value d(τ 0 )/N . Such a 3D distribution of x and k will be a Monte Carlo estimate of the function D(x, k, τ ) for a given value of τ with a statistical error proportional to 1/ √ N . Since a 3D distribution is difficult to visualise, in practice one usually makes 1D or 2D histograms of any combination of x and k. In addition, one can impose arbitrary cuts on any of these variables.
One can formally prove that the above algorithm gives a correct solution to Eq. (31), i.e. that the expectation value of a MC weight associated with a random walk trajectory, as described above, is equal to the function D(x, k, τ ). We skip such a proof here -it will be provided in our future publication dedicated to the MCMC algorithm and its implementation.
In the above MCMC algorithm we have assumed that all random variables can be generated according to the respective probability distribution functions using standard Monte Carlo techniques, preferably the analytical inversetransform method or its combination with the branching method. Among the integration variables in Eq. (31) the most problematic is the variable z because its probability distribution function ζ(z) is too complicated to be sampled with the above methods. Details on how to deal with this using a combination of the branching method with the importance sampling are given in Appendix A.
The MCMC algorithm for solving Eq. (21) is analogous to the above -one only needs to set w(q) = 0 and k n = 0, n = 0, 1, . . .

Differential method
Direct temporal numerical integration of the Eq. (7) is another approach we use. First of all, the spatial grid with the constant step-size x = 1/N is created to keep N grid-points with the solution of D(x, τ ) at each time-step τ . The generation of this grid seems to be quite involved in order to obtain reasonable numerical results and we shall discuss this issue later on. The integral on the right-hand side of the Eq. (7) is divided into two parts: the gain part and loss part Both of them are evaluated at every time step by a simple midpoint rule, i.e.
where z is the spatial step-size equal to x. Other advanced methods 1 did not yield significantly better results and this type of interpolation function is the fastest choice, which in turn allows us to use very dense grids (in fact, the grids with as many as N = 16,384 grid-points were used to obtain numerical results that are reasonably accurate and are presented in this paper). Due to the simple midpoint approximation we were able to keep computational time less than one day on a computer system with the i7 CPU. After the spatial approximation, we end up with the following system of ordinary differential equations: fourth-order accuracy and adaptive time-step regulator, see e.g. [38]. The time step correction is accomplished by the following rule: where δ tol = 10 −12 is the tolerance parameter used in the simulations and ε is the truncation error indicator computed in the last step of the algorithm. The initial condition used in the simulations is the analytical solution (8) at the very small time τ 0 = 10 −4 which we use as an approximation of the δ-function in the case of the simplified kernel. In the full kernel case, we use the following approximation of the δ-function: with ε = 6 · 10 −3 . The numerical solution is then advanced in time and reported. Spatial grids used during the simulations have to be very fine to obtain stable and meaningful results. The problem lies in the gain integral (32) where the arguments of K and D are reciprocal. When we use the equally distributed fixed grid and try to compute the gain integral with z from the finite subset of grid-points, we get arguments for D from unequally distributed grid points due to x/z, i.e. a bad approximation of the gain integral as most of the values will be taken from the region close to 0. A substitution does not help here as it will just switch the reciprocal values from one term to the other. To overcome this problematic behaviour, the very fine grid is needed that has an inevitable effect on computational times. One possible solution is the adaptive mesh refinement together with a smart distribution of the grid-points -this type of approach is still investigated.

Numerical results
We have implemented the MCMC algorithms described in Sect. 4 in the C-language program called MINCAS (the acronym for Medium-INduced CAScades) as two independent MC generators. First, we performed numerical tests of the algorithm for the solution of Eq. (21) by comparing it with the analytical formula of Eq. (8) and with the numerical differential method described in Sect. 5 for the case of the simplified z-kernel function, i.e. with f (z) = 1.
In our numerical calculations presented below we have used the following input parameters values: q min = 0.1 GeV, E = 100 GeV, n = 0.243 GeV 3 , The results for the evolution time values: t = 0.1, 1, 2 and 4 fm are presented in Fig. 1. We can see a very good agreement between the three solutions: by the MCMC algorithm of MINCAS, by the analytical formula (8) obtained in [39] and by the differential method described in Sect. 5. The resulting distributions feature the turbulent behaviour, i.e. the energy is transported from the large-x region to the low-x region without accumulating in the intermediate values of x.
In Fig. 2 we show similar results as above, but for the exact z-kernel function as given in Eq. (2). The agreement between MINCAS and the differential method is similar as in Fig. 1 which confirms that our numerical solutions of Eq. (7) are also correct for the exact z-kernel. Of course, now the analytical solution is away from both of them because it works only for the simplified z-kernel -it is shown only for reference. One can see that the x distribution for the exact z-kernel differs considerably from the one for the simplified z-kernel, particularly in the region of the intermediate x values -the turbulent behaviour of the exact solution is stronger than of the approximate one.
Then, we performed tests of the MCMC algorithm for the for the x and k evolution of Eq. (19) implemented in MINCAS.
Since the integration over k of Eq. (1) gives Eq. (7), our first test was to check if using the algorithm for the x and k evolution of Eq. (19) we can reproduce the x distributions generated by the simpler algorithm for x-only evolution of Eq. (21). For this purpose we have produced inclusive histograms of x, i.e. without any restrictions on k. The results in the case of the exact z-kernel function for the evolution time values: t = 0.1, 1, 2 and 4 fm are shown in Fig. 3. As one can see, they are in a perfect agreement. This is an important, non-trivial test of the MCMC algorithm for the x and k evolution of Eq. (19) and its implementation in MINCAS, showing that it produces the correct x distribution.
Unfortunately, we could not make comparisons of the k distributions with the differential method because it turned out to be inefficient in solving the general evolution Eq. (1). Therefore, in the following we present a few figures with the results from MINCAS only, to show how the the mediuminduced QCD evolution affects transverse gluon momenta.
In Fig. 4 the k x and k T = k 2 x + k 2 y distributions are shown. These results have been obtained with w(q) of Eq. (6) for the evolution time values: t = 0, 0.1, 1, 2 and 4 fm. One can observe fast broadening of a very narrow initial Gaussian distribution of k with the increasing evolution time as well as the departure from the Gaussianity of the subsequent distributions. This non-Gaussian behaviour of the transverse momen- tum distributions can be explained by inspecting Eq. (20). As one can see, k n is the sum of n + 1 random variables. From the Central Limit Theorem it follows that for a fixed value of n a distribution of the random variable k n would converge to the Gaussian distribution. However, in this case n is also a random variable as it corresponds to the length of the randomwalk trajectory in the MCMC algorithm described in Sect. 4. The final distribution of k results from summing of all such trajectories, therefore it is not a single Gaussian distribution but a sum of an arbitrary number of Gaussian distributions with the same mean values and different widths (variances). Generally, the longer trajectory results in the larger width as the transverse momentum broadening due the mediumcollisions seems to dominate, for a given form of the function w(q), over its shrinking due to the emission branchings. As the evolution time increases the trajectories get longer and more Gaussian distributions with larger widths contribute to the overall k distribution, making it wider -this we observe in Fig. 4.
In Fig. 5 we show the dependence of the k T mean value on the x variable for the q-kernel function w(q) of Eq. (5) (the LHS plot) and for that of Eq. (6) (the RHS plot). One can observe increasing k T in the course of the evolution in the whole x region for small evolution times and its accumulation in the low-x region for large evolution times. This pattern is very similar to the one presented in Ref. [33] for some approximate analytical solution. One should also comment on a distinctive feature of the slope of the final-state cascades k T as a function of x as compared with the behaviour of k T as a function of x in the initial-state cascades [40]. Namely, for the final states one can see that the lower x the typical k T is lower, while for the initial-state cascades the opposite happens. From Fig. 5 it can also be seen that k T rises faster with the evolution time for w(q) of Eq. (5) than for that of Eq. (6). This can interpreted as a more efficient quenching by the non-equilibrated plasma than by the equilibrated one, however the shapes are very similar, so the rate of the quenching is similar.
Finally, in Fig. 6 we show examples of 2D distributions of k x vs. k y (upper row) and x vs. k T (lower row) for the exact z kernel and w(q) of Eq. (6). The LHS plots present initial distributions, i.e. for t = 0, while the RHS ones the evolved distributions at t = 2 fm. One can observe how the initial gluon distributions get 'diffused' in x and k in the course of the medium-induced QCD evolution. The apparent departure from the Gaussian k distribution can be clearly seen in the upper-right plot. In the lower plots the turbulent behaviour of the distribution in the x direction, as discussed above, is also visible.

Summary and outlook
In this paper we have obtained numerical solutions of the equations describing the inclusive gluon distribution as produced by a jet the propagating in QGP, given in Ref. [31]. These equations were reformulated as the integral equations which allows for their efficient solution using the newly constructed Markov Chain Monte Carlo algorithms implemented in the dedicated Monte Carlo program MINCAS. The results for the energy density (the x distribution) were cross-checked with algorithm based on a direct numerical solution of the integro-differential equation by applying the Runge-Kuttabased method, and for the simplified emission kernel also with the exact analytical solution [31]. The MCMC method turns out to be far more efficient in solving the above equations than the differential method.
The resulting distributions of the gluon density as function of the transverse momenta show some new features, not studied so far in the literature on this subject, i.e. the departure, as the evolution time passes, from the initial Gaussian distribution. This is a result of the exact treatment of the gluon transverse-momentum broadening due to an arbitrary number of the collisions with the medium together with its shrinking due an arbitrary number of the emission branchings. We observe this behaviour for two different forms of the collision kernel w(q).
In the future, we plan to study in a more detailed and systematic way a relation of our MCMC solution to the existing approximate solutions as well as to test other possible forms of the collision kernel w(q) and the quenching parameterq resulting from them (in the present study, in order to have a correspondence to existing results, we have used the standard value ofq = 1 GeV 2 /fm). This will allow to see how universal the pattern of the gluon distribution in QGP is. For instance, one can use some AdS/CFT models to obtain w(q). One can also use our MCMC-based method to solve more general versions of Eq. (1) or an even more general kinetic equation (which assumes thermalisation of soft gluons) obtained in Ref. [8], and perform a full parton-shower simulation of the final state based on the generated distribution. Data Availability Statement This manuscript has no associated data or the data will not be deposited [Authors' comment: The paper is a theoretical study and no data is currently available.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Branching method with importance sampling
As was said in Sect. 4, the random variable z cannot be easily generated according to the pdf ζ(z), because it is a complicated function. For this purpose we can utilise the importance sampling technique, i.e. we can replace ζ(z) with some simplerζ (z):ζ and compensate for this simplification with an appropriate MC weight. The above simplification affects, however, generation of the random variable q i because z has the joint pdf with q, namely ξ(z, q) given in Eq. (25), and also generation of τ i . In order to describe this in detail, let us introduce some useful notation: where W is given in Eq. (12). Then, we can express the product of probability densities of the random variables τ i , z i and q i in terms of the above functions: wherẽ and v(τ i , z i ) is the compensating weight for simplifications done in generation of the random variables τ i and z i : with It turns out that for the difference of the integralsκ( ) and κ( ) we can take the limit → 0 -this limit is finite and can be computed (e.g. numerically) once for a given function f (z), independently of ; the above value corresponds to f (z) given in Eq. (2) (e.g. for a simple case of f (z) = 1: = 0). This suggests that in the MC generation we can avoid calculation of the complicated -dependent integral κ( ), instead we can replace it with the simple integralκ( ) and compensate for their difference with MC weight of Eq. (49).
One can prove that the expectation value of the above weight corresponds to the solution of Eq. (19), i.e. D(x, k, τ ) .
In actual MC computations, the above expectation value is estimated (according to the Law of Large Numbers) by the arithmetic mean of the event-weight values for a given MC sample. Its statistical error is proportional to 1/ √ N , where N is the number of generated MC events. We skip the proof of Eq. (56) here -it will be given in our future publication devoted to details of the MCMC algorithm.
In the case of the algorithm for solving Eq. (21), the above method simplifies to the pure importance sampling of z without branching into q, and k i = 0, i = 0, 1, . . .