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 the 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 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 operational the jet quenching can be observed at much higher available energies in collisions of lead nuclei [6]. Still one of the open problems is to understand the details of 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 weak-coupling regime of QCD [7][8][9][10][11][12][13][14][15], or AdS/CFT where one assumes the plasma to be strongly coupled [16,17] or classical field theory based approach [18] (for reviews see [7,[19][20][21][22]). Some of the mentioned formalisms are implemented in Monte Carlo event generators [23][24][25][26][27][28].
In this paper we look closer at results obtained in [29,30] 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 [29][30][31] 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 [30,32,33] where for instance part the equation leading to broadening of transversal 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 [32]. However, we find that the exact numerical solution is not Gaussian. Furthermore, the numerical method allows to test consequences of assumptions about properties of the medium for distribution patterns of the jets emitted from the hard jet.
The paper is organised as follows. In Section 2 we introduce and overview of the equations for the jet energy distribution and for the inclusive gluon distribution. In Section 3 we present reformulation of the above equations making use of Sudakov-like 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 Sudakov-type form factor. Then, we give formal iterative solutions of these equations. In Section 4 we propose a Markov Chain Monte Carlo (MCMC) algorithms for numerical solutions of the above equations. In Section 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 Section 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 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 Section 7. Finally, in Appendix A we provide some further details on the MCMC algorithms.

Evolution equations
The evolution equation for gluon transverse-momentum-dependent distribution D(x, k,t) in the dense medium, obtained under the assumption that the momentum transfer in the kernel is small, reads [30] 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 longitudinal momentum fraction, k = (k x , k y ) -its transverse-momentum 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 [30] w(q) = 16π 2 α 2 s N c n q 4 , 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 [34] w(q) = 16π 2 α 2 s N c n 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 [30] This integro-differential equation can be solved analytically for a simplified case of f (z) = 1 and D(x,t = 0) = δ(1 − x) [30]: 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 → 0 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 cut-off 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 be then transformed into the following 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 function D(x 0 , τ 0 ) and D(x 0 , k 0 , τ 0 ), respectively. In Eqs. (19) and (21) there is an 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: 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 sampling -we 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 [35]. 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 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 inverse-transform 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.

Differential method
Direct temporal numerical integration of the Eq. (7) is another approach we used. 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 = 16384 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: with x i = (i + 0.5)/N and z j = ( j + 0.5)/N, i, j ∈ {0, . . . , N − 1}. The RHS in Eq. (35) is then used to advance the solution in time by the five-stage Runge-Kutta-Merson method with fourth-order accuracy and adaptive time-step regulator, see e.g. [37].
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. 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 behavior, 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 Section 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 Section 5 for the case of the simplified z-kernel function, i.e. for f (z) = 1. The results for the evolution time values: τ = 0.1, 0.2, 0.5 and 1.0 are presented in Fig. 1. We can see a perfect agreement between the results from MINCAS and the analytical solution as obtained in [38]. The results from the differential method agree very well with the other ones for mildly varying parts of the x distribution, while in the parts with strong slopes they are slightly shifted to the right -this is due to limitations of this method, as explained in Section 5. But the overall agreement is satisfactory given the complex form of the equation. The resulting distributions feature the turbulent behaviour i.e. the energy is transported from the large x to the low x 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. 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 observe that the x distribution for the exact z-kernel departs more and more from the solution for the simplified kernel with the increasing evolution time.
Then, we performed tests of the MCMC algorithm for the for the x and k evolution of Eq. (19) implemented in MINCAS using the following input parameters values: E = 100 GeV, n = 0.243 GeV 3 ,q = 1 GeV 2 /fm .
Since 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 equation (1). Therefore, in the following we present a few figures with the results from MINCAS only, to show how the the medium-induced 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.
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) with q min = 0.1 GeV (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. [32] 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 [39]. 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 also 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 mediuminduced 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. These equations were reformulated as the integral equations which allows for 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-Kutta method, and for the simplified emission kernel also with the exact analytical solution. 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 to the collisions with the medium together with its shrinking due the emission branching. 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 quenching parameters resulting from themq (in the present study, in order to have a correspondence to existing results, we have used the standard valuê q = 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 the equation (1) or perform a full parton-shower simulation of the final state based on the generated distribution. 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: and v(τ i , z i ) is the compensating weight for simplifications done in generation of the random variables τ i and z i : with where ∆ = lim ε→0 [κ(ε) − κ(ε)] ≈ 3.57066164.
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. (48). Because of the change in the τ-variable pdf: ρ(τ) →ρ(τ), we also need to modify accordingly the stopping rule: If the initial density D(x 0 , k 0 , τ 0 ) is a complicated function, we can approximate it with some simpler functionD(x 0 , k 0 , τ 0 ), construct the corresponding pdf: and apply the compensating weight u(x 0 , k 0 ) = D(x 0 , k 0 , τ 0 ) D(x 0 , k 0 , τ 0 ) .
One can prove that the expectation value of the above weight corresponds to the solution of Eq. (19), i.e. E[w γ (x, k, τ)] = 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. (55) 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. . . ..