Medium induced QCD cascades: broadening and rescattering during branching

We study evolution equations describing jet propagation through quark-gluon plasma (QGP). In particular we investigate the contribution of momentum transfer during branching and find that such a contribution is sizeable. Furthermore, we study various approximations, such as the Gaussian approximation and the diffusive approximation to the jet-broadening term. We notice that in order to reproduce the BDIM equation (without the momentum transfer in the branching) the diffusive approximation requires a very large value of the jet-quenching parameter q̂\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{q} $$\end{document}.


Introduction
Quantum Chromodynamics (QCD) is the very well established theory of strong interactions with rich structure and many phases [1]. Here we want to focus on a jet-quenching phenomenon predicted in [2,3], and observed experimentally at RHIC [4] and LHC [5]. The jet quenching is a suppression of propagation of jets in quark-gluon plasma (QGP) due to jet-plasma interactions. This process has many phases, recently discussed in refs. [6,7], see also [8]. The jet-quenching phenomenon is approached from many directions: the kinetic theory [9][10][11][12][13][14][15][16][17][18][19], Monte Carlo methods [20][21][22][23][24][25], the AdS/CFT [26]. Furthermore, it is a multi-scale problem which, however, allows for factorisation in time. In particular, according to refs. [6,7], in the first phase the jet propagates according to the vacuum-like parton shower with ordering in an angle, while in the next stage the coherence is broken and jet propagates through plasma experiencing elastic scatterings and branching -in this stage there are many soft radiations and wide-angle emissions. In the last stage, when jet leaves medium, again the vacuum-like emissions dictate its time evolution. In this paper, we focus on the second phase of the jet propagation through QGP In particular, we investigate what is the contribution of momentum transfer during branching to the broadening pattern of the jet. To address this problem, we solve the equation proposed in [27,28] which is a generalised version of the equation solved by three of us in ref. [29]. 1 In this approach, QGP is modelled by static centres and the jet interacts with it weakly, jet propagating through plasma branches according to BDMPS-Z mechanism [9][10][11][12][13][14][15][16][17][18] and gets broader due to elastic scattering with plasma.
In ref. [29] it has been observed that accounting for the broadening term beyond the diffusive approximation leads to the non-Gaussian broadening for jet observables. It turns JHEP04(2021)014 out that the non-Gaussianity leads to much stronger broadening of the cross section for decorrelations than the Gaussian approximation.
In the current study, we propose detail study of impact of momentum exchange during branching and its contribution to the broadening. Experimentally, the broadening is rather small and its observation at the LHC energies is hindered by the vacuum effects [31,32]. There are also possible effects which could give negative contribution to the broadening [33][34][35]. Furthermore, a more realistic model of the medium accounting for its expansion [34,36,37] will probably reduce the amount of the broadening. While at present we do not account for a more realistic scenario, i.e. the expansion of the medium we mimic it by scaling theq parameter. We observe that reducing its value leads to the smaller broadening.
The paper is organised as follows: in section 2 we discuss the version of the transversemomentum-dependent BDIM (Blaizot, Dominguez, Iancu, Mehtar-Tani) equation [28] where the momentum transfer in the kernel is taken into account and we present its solution with the use of Monte Carlo methods. In section 3 we compare the BDIM equation to some of its approximations, i.e. the case where transverse momentum in the branching kernel is neglected, the case when the broadening term is represented by the diffusive approximation, and the Gaussian approximation where the transverse momentum and the longitudinal momentum are factorised. We conclude our work in section 4. In appendix A we present one of the Monte Carlo algorithms for solving the full BDIM equation, 2 while in appendix B we describe a numerical method used to solve the diffusive approximation of the BDIM equation.

Momentum-transfer-dependent BDIM equation and its solution
The evolution equation for the gluon transverse-momentum-dependent distribution D(x, k, t) in the dense medium reads [28] (2.1) The kernel K(Q, z, xp + 0 ) which accounts for the momentum-dependent medium induced branching is given by The other one is an extension of the algorithm employed in the Monte Carlo program MINCAS, described in ref. [29], and will be presented elsewhere.

JHEP04(2021)014
. (2.4) where p + 0 ≡ E is energy of jet entering the medium, x -is longitudinal momentum fraction of mini jet, k = (k x , k y ) -is transverse-momentum vector of mini jet,q -the quenching parameter, α s -the QCD coupling constant and N c -the number of colours.
The elastic collision kernel C(l) is given by where the function w(l) models the momentum distribution of medium quasi-particles. We consider two scenarios: 1. The out-of-equilibrium distribution [28]: with l = (l x , l y ) being transverse-momentum vector and n -the density of scatterers.
2. The situation where the medium equilibrates and the transverse-momentum distribution assumes the form obtained from the Hard Thermal Loops (HTL) calculation [38].
In this case the medium is characterised by a mass scale given by the Debye mass m D : , The equation (2.1) has been solved using the Monte Carlo program MINCAS by extending the algorithm presented in [29] (to be described elsewhere) and, independently, using another Monte Carlo algorithm described in the appendix A. The two solutions have been checked to be in a good numerical agreement. Here we present the results from MINCAS obtained using the following input parameters: In figure 1 we show the k T distributions as well as k T as a function of x for the evolution time values t = 0, 0.1, 1, 2, 4 fm. The detailed discussion of the solution is presented in the next section where we also discuss comparisons to the approximations of the BDIM equation.   Figure 1. The k T and k T vs. log 10 x distributions for the evolution time values t = 0, 0.1, 1, 2, 4 fm, for the full emission kernel K(Q, z, p + ) (denoted as K(Q, z)) and the collision term of eq. (2.6).

Comparison of BDIM equation to its approximations
In this section we will discuss various approximations of the BDIM equation.
• The first approximation that we consider is the case when the momentum transfer during branching is neglected. In this case, as demonstrated in ref. [27], the branching kernel simplifies to a purely collinear one and the transverse momentum dependence comes basically from the elastic scattering. The equation reads • One can further simplify the BDIM equation by expanding the elastic collision term and using the diffusive approximation [27] to obtain In the above equation, as compared to ref. [27], we have neglected the mild logarithmic dependence ofq in the diffusion term on k T .
• Eq. (3.3) was also solved approximately in ref. [39]. To arrive at the solution, the branching term was neglected and the Gaussian ansatz was used. The solution reads where In the above, it is assumed that k 2 ⊥ < ω 2 = (xE) 2 , and the parameters are:ᾱ s = 0.3, Similarly to eq. (2.1), eq. (3.1) has been solved using the Monte Carlo programs, basically re-obtaining the result from ref. [29], while eq. (3.3) has been solved with the help of the numerical method described in appendix B. For eq. (2.1) and eq. (3.1) we have used both the functions w(l) from eqs. (2.6) and (2.7) to describe the collision term when using both the full kernel and the simplified one. In the case of the full kernel, we have also performed calculations without the collision term. For all the presented results, we have used the parameters given in the previous section.
In figure 2 we show the k T distributions of the six cases studied for evolution time values t = 0, 0.1, 1, 2, 4 fm. We directly notice that the Gaussian approximation fails to describe any of the other results. The nearest distribution is the one with the full kernel and no collision term, which approaches a Gaussian shape, but with a much wider width. The other distributions (with the collision term) show fast broadening of the initial Diracδ-like distribution, exhibiting the non-Gaussian shape. The broadening is faster with w(l) given by eq. (2.6) than with the one given by eq. (2.7), i.e the broadening is faster with out-of-equilibrium momentum distributions of the medium quasi-particles.
In figure 3 we present the dependence of the mean value of k T on log 10 x. For all cases, k T grows with time and with x. It is still true for the Gaussian approximation, even if the distribution for the different evolution time join each other under certain values of x. We can clearly see in these figures a different behaviour around x = 1 between the JHEP04(2021)014 Figure 3. The k T vs. log 10 x distributions for the evolution time values t = 0, 0.1, 1, 2, 4 fm, for different kernels: the Gaussian approximation, K(z) and K(Q, z, p + ) (denoted as K(Q, z)), and different collision terms: no collision term, the collision term as in eq. (2.6) and in eq. (2.7), respectively.
distributions corresponding to the z-only dependent kernel and the ones corresponding to the full kernel which show a drop. This drop results from the fact that the evolution starts at x = 1 with k T = 0 and already a single soft emission with the Q-dependent kernel K gives to the emitter a significant k T -kick, which is not the case for the z-only dependent emission kernel. This effect is more pronounced for the shortest evolution time t = 0.1 fm, because in this case the (x, k T )-distribution is strongly peaked at x = 1 and k T = 0, while for the longer evolution times this peak is smeared out, so the contribution from x = 1 and k T = 0 to k T is much smaller. Except for the drop near x = 1 for short evolution times with the full emission kernel, the k T distributions increase with x. The Gaussian approximation gives the lowest k T values, while they are the highest for the evolution with the full emission kernel -and these, in particular, are higher than in the case with the z-only dependent emission kernel. This results from the fact that in the former case the k T -broadening is produced not only in the collisions with the medium (due to the C(l) term in eq. (2.1)), but also in the emission process (due to the Q-dependence of the kernel K in eq. (2.1)).
In figure 4 we present distributions integrated over the transverse momenta for four values of the evolution time: t = 0.1, 1, 2, 4 fm. We see that all the transverse-momentumdependent distributions, as a consequence of momentum conservation, collapse to the same x-dependent distributions. This further confirms that the study of the transversemomentum dependence allows for more detailed study of the dynamics of the branching process. It also constitutes an important numerical cross-check that all our algorithms for the transverse-momentum-dependent evolution satisfy the condition: D(x, k, t) .  An interesting question is what is the domain of applicability of the diffusive approximation that was used in order to reduce eq. (3.1) to eq. (3.3). The approximation is advocated as a systematic expansion around k T that should be valid for rather low values of k T . However, from the explicit solution in figure 5 we see that the solution of the eq. (3.3) is reasonably reproduced in the diffusive approximation if we allowq to be very large. This actually is in agreement with the interpretation ofq as the average transverse momentum. Therefore, we conclude that one can describe large transverse momentum using just the diffusion approximation, but one should allow this new effectiveq to be large and different from the one in the complete equation. We also see that the diffusive approximation with the standardq preserves the general pattern of eq. (3.1), but is much narrower than the solution of the equation before the expansion. This feature is better visible in the plot of the k T as a function of x which we show for different values ofq as well as for different values of t. From these results we conclude that, while the diffusive approximation is qualitatively fine, it is rather crude quantitatively.
To complete the analysis of the k T spectrum we study on figure 6 its dependence on q for the three cases ofq = 0.5, 1, 2 GeV 2 /fm. We see that, in general, it is not a trivial dependence, in a sense that increasingq will just broaden the distribution. This is the case only for the Gaussian approximation and w(l) = 0. The interpretation of this is the JHEP04(2021)014 In the diffusive approximation σ k0 = 1 GeV was used; note also that the evolution times forq = 1500 GeV 2 /fm are equal to t = 0, 1, 2, 3 fm (τ ≡ t/t * = 0.0675 when t = 1 fm) in the case ofq = 1 GeV 2 /fm.

JHEP04(2021)014
following. In these two casesq enters to some extent trivially: in the former case as a factor modifying k 2 br , while in the latter in the branching term only. In the remaining cases,q which controls the broadening enters the branching kernel and is hidden in both the branching term and the elastic scattering term. Interplay of these two effects results in the structure visible in these cases.

Conclusions and outlook
We have solved and studied the BDIM equation as well as its various approximations, i.e. the no-momentum-transfer approximation, the diffusive approximation and the Gaussian approximation. We conclude that the momentum transfer during branching gives additional broadening that is non-negligible. Furthermore, the diffusive approximation of the elastic scattering kernel is a rather crude approximation to the BDIM equation. In the future it will be interesting to investigate the case of the expanding medium as well as to account for coupled evolution of quarks and gluons. Furthermore it will be interesting to see the signature of the rescattering during branching in some final state. One of the possibilities is to study decorelations of jets following [40].

A Monte-Carlo algorithm
With the help of the Sudakov form-factor where the notation |q| > q ↓ should indicate that the integration runs over all q except those where |q| < q ↓ , it can be shown that the following integral equation is equivalent to the integro-differential equation eq. (2.1): in the simultaneous limits of → 0 and q ↓ → 0.

JHEP04(2021)014
The individual terms in eq. (A.2) can be associated with probabilities: • The probability that the fragmentation function at time t gets a contribution from the fragmentation function at time t without additional splitting or scattering between t and t (but at t and t some particle interaction occurs): • The probability density that the fragmentation function at the momentum fraction x and the transverse momentum k gets a contribution from the fragmentation function at the earlier time t at the momentum fraction y and the transverse momentum q via a splitting with momentum fraction z and transverse momentum Q, where x = zy and k = Q + zq: Thus, the probability for a splitting with a certain z value (independent of the value of Q) is given as zK(z) • The probability density that the fragmentation function at the transverse momentum k gets a contribution from the fragmentation function at the earlier time t at the transverse momentum q via a scattering with the transverse momentum Q, where k = Q + q: .
Thus, it is possible to obtain solutions for eq. (2.1) via a Monte-Carlo algorithm, where a distribution D(x, k, t) that obeys eq. (A.2) can be obtained by selecting independently of one another a large number N ev of sets (x, k), which follow D(x, k, t).
In each of the N ev cases, the x and k values are obtained in the following way: • Some initial values x 0 , k 0 are set together with the time t 0 of the start of the evolution.
• For every set ( • The previous step is repeated until for some time t j j ∈ N, it is found that t j ≥ t. Then the algorithm gives x = x j−1 , k j = k j−1 and stops.

JHEP04(2021)014
The selection of a set (x i+1 , k i+1 , t i+1 ) from a set (x i , k i , t i ) is done in the following way: 1. Select time t i+1 of next splitting/scattering by first choosing a random number R ∈ [0, 1] from a uniform distribution and then solving the equation .
The result of this calculation is 2. Determine whether a splitting or scattering occurs. This is done, by first selecting a random number R ∈ [0, 1] from a uniform distribution. If a scattering occurs, otherwise a splitting.
3. If a splitting occurs, determine x i+1 and k i+1 as follows: (a) Select z from K(z) by choosing a random number R ∈ [0, 1] from a uniform distribution and then solve the equation This equation is solved approximately by first tabulating values of z 0 dz z K(z ) for a set of z values that is sufficiently dense for the desired accuracy and then searching from this table the z value, which is the closest to the one that solves eq. (A.11). After selection of a, the value of Q = 2k 2 br a is calculated. While the values of a can assume any positive value, we here constrain the values to the region a ∈ [0, π] in order to avoid the region where sin(a)e −a becomes negative. Indeed the splitting function in the form of eq. (2.2) was deduced in ref. [27] in the harmonic approximation, which needs corrections at large momentum scales.

JHEP04(2021)014
The terms on r.h.s. of the eq. (B.1) are evaluated by central differences (the Laplacian of k with one-sided approximations at the boundaries of the computational domain) and by the box-rule (the integral term): (B.4) A numerical grid is equidistant and 2-dimensional (we drop the φ-dependence due to the symmetry of the problem): We solve eq. (B.4) to obtain the functions D i,j (t n ) = D(x i , k j , t n ) at given points x i , k j and a time level t n . The initial condition is given by eq. (B.2). The number of grid points for x and k is increased up to N x = 10240 and N k = 1000 with x ∈ [0, 1] and k ∈ [0, 50] (k max = 50) for the case ofq = 1500 GeV 2 /fm, for otherq we used coarse grid with N x = 1024 and N k = 200.
We use a fourth-order Runge-Kutta method to obtain the numerical solution of the eq. (B.4) in time (the Cash-Karp method with the adaptive time stepping [41] is employed). The time step is being changed according to the following formula: where TOL = 10 −6 is a tolerance and E is the maximal error in the last step of the embedded Runge-Kutta method. In order to minimise the computational time, the numerical code was parallelized and implemented in NVIDIA CUDA (double precision was used in computations).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.