Anomalous diffusion in QCD matter

We study the effects of quantum corrections on transverse momentum broadening of a fast parton passing through dense QCD matter. We show that, at leading logarithmic accuracy the broadening distribution tends at late times or equivalently for large system sizes $L$ to a universal distribution that only depends on a single scaling variable $k^2_\perp/Q^2_s$ where the typical transverse momentum scale increases with time as $\ln Q_s^2 \simeq (1+2 \beta ) \ln L - \frac{3}{2}(1+\beta )\,\ln\ln L$ up to non-universal terms, with an anomalous dimension $\beta \sim \sqrt{\alpha_s} $. This property is analogous to geometric scaling of gluon distributions in the saturation regime and traveling waves solutions to reaction-diffusion processes. We note that since $\beta>0$ the process is super-diffusive, which is also reflected at large transverse momentum where the scaling distribution exhibits a heavy tail $k_\perp^{4-2\beta }$ akin to L\'{e}vy random walks.

We study the effects of quantum corrections on transverse momentum broadening of a fast parton passing through dense QCD matter. We show that, at leading logarithmic accuracy the broadening distribution tends at late times or equivalently for large system sizes L to a universal distribution that only depends on a single scaling variable k 2 ⊥ /Q 2 s where the typical transverse momentum scale increases with time as ln Q 2 s (1 + 2β) ln L − 3 2 (1 + β) ln ln L up to non-universal terms, with an anomalous dimension β ∼ √ αs. This property is analogous to geometric scaling of gluon distributions in the saturation regime and traveling waves solutions to reaction-diffusion processes. We note that since β > 0 the process is super-diffusive, which is also reflected at large transverse momentum where the scaling distribution exhibits a heavy tail k 4−2β ⊥ akin to Lévy random walks.
Transverse momentum broadening (TMB) of energetic quarks and gluons while traversing QCD matter plays a central role in a variety of processes studied at colliders to probe QCD, ranging from jet suppression in heavy ion collisions [1, 2] to transverse momentum dependent gluon distribution function that encode information on the 3D structure of the proton and nuclei in high energy collisions in particular at small Bjorken x where gluon saturation is expected to take place [3][4][5].
High energy partons experience random kicks in hot or cold nuclear matter that cause their transverse momentum k ⊥ w.r.t. to their direction of motion to increase over time. The leading order elastic process is given by a single gluon exchange via Coulomb scattering and leads to an approximate brownian motion in transverse momentum space, so long as the mean free path is larger than the in-medium correlation length, where the typical transverse momentum square scales linearly with system size L, namely k 2 ⊥ typ ∼q L, whereq is the diffusion coefficient.
Moreover, radiative processes can also increase transverse momentum of the leading parton due to recoil effects. This question was recently addressed, and it has been shown that such contributions, albeit suppressed by the coupling constant α s , are enhanced by double logarithms which must be resummed to all orders when α s ln 2 L ∼ 1 [6][7][8][9].
In this letter we go beyond this result by investigating in more detail the consequences of the non-local nature of the quantum corrections on the TMB distribution. We find in particular that the latter reaches a universal scaling solution at late times (large L) that we compute analytically along with its sub-asymptotic deviations exploiting a formal analogy between the present problem and traveling wave solutions to reaction-diffusion processes [10]. As a consequence of the self-similarity characterizing the anomalous random walk, the TMB distribution is of Lévy type. It is in particular associated with a heavy power law tail describing long rare steps which extends over a large range of transverse momenta above the typical transverse momentum scale.
Lévy flights are ubiquitous in nature and span a wide variety of stochastic processes in biological systems [11,12], molecular chemistry [13], optical lattice [14], turbulent diffusion and polymer transport theory [15,16]. Furthemore, heavy tailed distributions are also observed in self-organized critical states [17,18]. In this work, we point out for the first time another occurrence of such random walks in the transport of eikonal partons in dense QCD matter and we compute the anomalous exponents that characterize the deviation to standard diffusion.

QUANTUM CORRECTIONS TO TRANSVERSE MOMENTUM BROADENING IN QCD MEDIA
The TMB distribution is related to the forward scattering amplitude S(x ⊥ ) of an effective dipole in color representation R = A, F with transverse size x ⊥ (see e.g. [7,19,20]) via a Fourier transform, (1) Considering the dipole formulation in position space allows for a straightforward resummation of multiple scattering by exponentiating the single scattering crosssection, so long as the interactions between the dipole and the medium are local and instantaneous. Thus, we may write S( [21]. The latter relation defines the quenching parameter in the adjoint representation which is assumed to be a slowly varying function of x ⊥ . At tree-level it readŝ up to powers of x T µ suppressed terms. For a weakly coupled QGP the bare quenching parameterq 0 and the arXiv:2109.12041v1 [hep-ph] 24 Sep 2021 infrared transverse scale µ 2 [22] are related to the Debye screening mass in the QGP or to the inverse nucleon size in a nucleus. It is customary to define the emergent saturation scale Q s (L) via the relation S(x 2 ⊥ = 1/Q 2 s (L)) ≡ e −1/4 , or equivalently,q(Q 2 s (L))L ≡ Q 2 s (L). This definition is standard in small-x physics [23,24], and is also motivated by Molière theory of multiple scattering [25][26][27] in which Q s is the transverse scale that controls the transition between the multiple soft scattering and the single hard scattering regimes. Given Eq. (2), one finds Q 2 s ∼q 0 L ln(q 0 L/µ 2 ) at tree level. Beyond leading order in α s , one has to account for real and virtual gluon fluctuations in the effective dipole with lifetime τ smaller than the system size. Such fluctuations yield potentially large contributions of the form q (1) ∼ α sq (0) ln 2 (L/τ 0 ) with τ 0 L a microscopic scale of order of the mean free path [6]. These radiative corrections to the quenching parameter can be resummed to double logarithmic accuracy (DLA) via an evolution equation ordered in τ [6,8,9]: (3) enforces the gluon fluctuations to be triggered by a single scattering with plasma constituents whose contribution is logarithmically enhanced compared to multiple scattering. The final value of τ is fixed at the largest time allowed by the saturation condition, i.e. at the time scale τ such that Q 2 s (τ ) = k 2 ⊥ , as long as τ < L and τ = L otherwise. The latter case corresponds to the dilute regime [28], that is when k 2 ⊥ Q 2 s (L). In this letter, we address both analytically and numerically the non-linear system (3)-(4) [29]. Analytic solutions are in general difficult to obtain, however, a solution for the linearized problem that consists in approximating Q 2 s (τ ) q 0 τ for the emission phase space can be found [30,31]. For a constant initial condition it readŝ . Formally, this linearization is valid in DLA, which captures all the terms of the formᾱ n s Y 2n assuming Y ∼ ρ, but it misses sub-leading corrections of the formᾱ n s Y 2n−1 ln Y which are parametrically larger than the single logarithmic ones. This is one of the novelty of the present study, enabling us to highlight the geometric scaling property of the transverse momentum diffusion in QCD and to compute the scaling deviations.

GEOMETRIC SCALING AND TRAVELING WAVES
The TMB distribution is said to obey geometric scaling if it is only a function of k 2 ⊥ /Q 2 s (L) as a result of scaling invariance of the radiative process at late times. More precisely, we would have for some scaling function f to be determined. Note that the linearized version of Eq. (3) satisfies a similar scaling relation with the notable difference that the argument of The non-linearity of Eq. (3) enforces the evolution to be controlled by a single momentum scale Q s (L) [33].
Geometric scaling was extensively studied in the context of deep inelastic scattering, where it has been shown that the gluon distribution g(x, Q 2 ) at small x satisfies this property over a broad region of photon virtuality −Q 2 [34][35][36]. We shall demonstrate that TMB exhibits similar properties.
Traveling waves (TW) solution. Remarkably, it is possible to find the scaling function f for the non-linear problem defined by Eqs. (3)-(4). In terms of the variables Y and ρ, the differential equation satisfied byq reads with ρ s = ln(Q 2 s /(q 0 τ 0 )). Let us now look for a scaling solution of the formq Plugging this expression in Eq. (7), one gets the following differential equation In order for f to be a function of x only at large Y , the derivative dρ s /dY must converge towards a constant c, which can be interpreted as the speed of a traveling wave that propagates to the right on the ρ axis. This is reminiscent of the TW solutions [10] to the Balitsky-Kovchekov (BK) equation [37,38].
The initial conditions to this second order linear differential equation follow from the non-linear treatment of the saturation boundary: on the isoline x = 0, meaning ρ = ρ s (Y ), one has f (0) = 0 and f (0) = (c − 1)/c. It is then straightforward to solve Eq. (8) with the scaling form f (x) ∝ e βx where the slope β is solution of the quadratic equation −cβ 2 + (c − 1)β −ᾱ s = 0. For any initial condition such thatq(Y = 0, ρ) < e βρ when ρ → ∞ [39], which is the case in the present letter, the minimum value for β, that satisfies the additional constraint −2cβ + (c − 1) = 0, will be realized. This fixes the speed to c = 1 + 2 ᾱ s +ᾱ 2 s + 2ᾱ s 1 + 2 √ᾱ s . Hence, 10 −12 10 −9 10 −6 10 −3 10 0 for this critical value that minimizes the TW speed, the solution to (8) takes the form In terms of the physical variables the k 2 ⊥ dependence ofq that enters the broadening distribution reads, in the large L limit, which is continuous and derivable everywhere.
To make the interpretation of these results in terms of TW's more transparent we inserted Eq. (10) in S(x ⊥ ) and plotted the result in Fig. 1 for several values of L. We see that the traveling wave propagates from right to left (from large to small x T ) with increasing L. However, once plotted in terms of x 2 T Q 2 s as shown in the inset of Fig. 1, they all lie approximately on the same universal scaling function given by Eq. (1)-(10) (dotted black curve). The observed deviations will be discussed in what follows.
Sub-asymptotic corrections. We turn now to the calculation of the sub-asymptotic corrections to the geometric scaling solution (10). Near the wave front, typically for x 1, we can look for a solution of the form inspired by the traveling waves ansatz that solves the BK equation [40] and more generally FKPP-like equations [41,42]. Plugging this ansatz inside Eq. (7), one gets a differential equation for G(z). Because the coefficient of Y and Y α (assuming α > 0) in this equation must vanish, we recover the two previous constraints that fix the value of c and β. Then, neglecting the power suppressed terms Y −1 and Y −α−1 , one finds The homogeneity condition implies that the coefficient α must be equal to 1/2 so that the deviation from the scaling form near the wave-front grows in a diffusive way as Y increases.
For initial conditions such that β < β, the large z behavior of the function G(z) constrains the acceptable values of the coefficient in front of G(z) in this differential equation. In fact, as shown in [39,40], one must have This yields the solution The value of b we extract from this analysis is novel and a consequence of the non-linearity of the saturation boundary. In the linearized problem with the lower bound in the integral of Eq. (7)  scaling limit with its deviation provided by the function with ρ s (Y ) = cY + b ln(Y ) + const. This solution is independent of the initial condition (for physically relevant ones), and only depends on the value ofᾱ s via the coefficients c, β and b. The resummed TMB distribution displays a universal behavior independent of the non-perturbative modeling of the treelevel distribution often used as an initial condition for non-linear small x evolution [44,45]. It can therefore provide a model-independent functional form for the initial condition of the BK equation, that includes gluon fluctuations enhanced by double logs,ᾱ s ln 2 A 1/3 , inside the nucleus target to all orders.

SUPER-DIFFUSION AND MODIFICATION OF RUTHERFORD SCATTERING
In this last section, we investigate the physical consequences of the scaling solution (10) forq(k 2 ⊥ ) on the TMB distribution given by Eq. (1), in particular at large k T , where the distribution is characterized by rare events that are sensitive to the point-like nature of the medium scattering centers [46]. First, it is straightforward to see that in the large L limit, the TMB distri-bution P(k T ) is only a function of k T /Q s (L). In Fig. 2, we show the distribution as a function of k T /Q s (L) with Y = ln(L/τ 0 ) = 4, for the following set-ups: (i) in dashdotted grey, at tree-level, using Eq. (2), (ii) after quantum evolution obtain by numerically solving Eqs. (4) in red, (iii) in blue, using the expression Eq. (16) that includes sub-asymptotic corrections to the scaling limit, (iv) finally, in dashed black, the scaling limit Y → ∞ of Eq. (16). Interestingly, the sub-asymptotic corrections account for the relatively large deviations between the asymptotic curve and the exact numerical result at the moderate value of L = 6 fm.
The k T distribution exhibits two different regimes: the region of the peak, near Q s (L) and the large k T tail, with k T Q s (L). These results can be interpreted in term of a special kind of random walk (here in momentum space) called Lévy flight. Such a remarkable connection with statistical physics enables us to highlight some interesting features (i) self-similar dynamics (ii) super diffusion (iii) power-law tail with slower decay than the Rutherford k −4 ⊥ behavior. In order to further the connection with the physics of anomalous diffusion, consider the scaling limit of the TMB distribution in the vicinity of the peak where the shape of the distribution is controlled by the first line in Eq. (10). Using this solution, one finds that S(x ⊥ ) In momentum space, it implies that the distribution P(k ⊥ , L) satisfies a generalized Fokker-Planck equation, ∂P/∂L ∝ −(−∆) 1−2β P, where the so-called fractional Laplace operator (−∆) γ/2 is defined by its Fourier transform |x ⊥ | γ [47,48]. This fractional diffusion equation (without external potential) is satisfied by the probability density for the position of a particle undergoing a Lévy flight process in two dimension [49] with stability index γ = 2 − 4β 2 − 4 √ α s + O(α s ).
Because of its heavy tail (to be discussed thereafter), the mean k 2 T of the TMB distribution is not defined. Nevertheless, it is possible to introduce a measure of the characteristic width of the k T distribution, and study its behavior as a function of the medium size L. In what follows, we shall use the median value k T med of k T P(k T ) [50] which is shown in Fig. 3 for three different scenarios. In grey, we plot the tree-level resulting from Eq. (1) and (2). The median scales approximately like (L ln L) 1/2 , which up to the logarithmic factor resulting from the Coulomb logarithm in Eq. (2), exhibits the standard diffusion scaling. The red line is the median of the k T distribution obtained using the resummed value ofq with fixed coupling, after numerical resolution of Eqs. (4). We then compare this result with our analytic prediction (12) (assuming k T med ∝ Q s , which is represented in blue in Fig. 3 . Remarkably, the agreement is excellent down to rather small values of L ∼ 3 fm. Since c/2 > 1/2, the median grows faster than √ L at large L, illustrating the super-diffusive behavior of TMB beyond leading order, with a deviation to the standard diffusion of order √ᾱ s . Heavy-tailed distribution. Another important aspect of Lévy flights is the power law decay of the step length distribution for a Lévy walker [52,53]. This reflects the fact that long jumps with arbitrary length may occur with non-negligible probability. In the problem at hand, this power-law tail can also be understood as a consequence of the self-similar nature of overlapping successive gluon fluctuations. The tail of the TMB distribution is controlled by the large k 2 ⊥ behavior ofq(k 2 ⊥ ), and consequently, by the exponential in the second line in (10). Note, however, that the scaling limit encompasses two stability indices: one controlling the peak and the median, as discussed above, while the other controls the tail of the distribution (cf. (10)). Without loss of generality, one can derive the leading behavior of P(k ⊥ ) at large k ⊥ by expanding the dipole S-matrix for small dipole sizes, as a result the Fourier transform can be approximated by [54]: up to logarithmically suppressed terms. This formula quantifies the deviations from the Rutherford scattering cross-section that are induced by radiative corrections. Applying Eq. (17) to our scaling solution (10), one finds the tail with ν = 2 − β + O(ln(x)/x) and x = ln(k 2 ⊥ /Q 2 s (L)).
The corrections to the power-law behavior are due to the pre-factor in the second line of (10). The power of the tail deviates from the tree level Rutherford behavior by ∼ −2 √ᾱ s . The form of ν is correct in the strict scaling limit L → ∞. For finite L values, the 1/k 4 ⊥ tail is recovered at very large k T , as can be inferred from Eq. (5) which yields ν = 2 − 2 ᾱ s Y /x (when x Y ). The fact that geometric scaling extends in the tail region is known in the context of saturation physics as the "extended geometric scaling window" corresponding to Q s k T Q 2 s /µ. Finally, note that this analysis is valid so long as E/k 2 T L, where E is the energy of the fast parton. This allows us to neglect the quantum diffusion of the dipole in the medium. In the opposite case the quantum phase is suppressed when E/k 2 T τ L and the evolution is expected to be of DGLAP type with the substitution ln L/τ 0 → ln E/(k 2 T τ 0 ) [55]. Conclusion. In summary, we have studied the transverse momentum distribution of a high-energy parton propagating through a dense QCD medium, including resummation of radiative corrections within a modified double logarithmic approximation which accounts for the non-linear dynamics due to multiple scatterings that restrict the phase space for quantum fluctuations. We have found that the non-linearity and self-similarity of overlapping multiple gluon radiations lead to a universal scaling limit at large system sizes, which exhibits a super-diffusive regime and a power-law decay akin to Lévy flights. Although at very high k ⊥ , the distribution is characterized by point-like interactions of Rutherford type, for moderately large k ⊥ we observe a weaker power due to the non-local nature of the interactions which is the hallmark of scale invariant phenomena.
Concerning phenomenological applications, we point out the relevance of our analytic solutions in the study of nuclear structure at high energy as it provides a new initial condition for non-linear evolution of the gluon distribution. We leave for future work the question of the experimental detection of this emergent QCD phenomenon in heavy ion collisions as well as running coupling corrections which are expected to yield mild scaling violations.

Numerical implementation of the evolution equation
Differential problem. We discuss here the numerical method used to solve the evolution equation forq with exact treatment of the saturation boundary for the gluon emission phase space. The evolution equation in differential form reads The existence of a solution is not guaranteed since the boundary of the k ⊥ integral depends onq(τ, k 2 ⊥ ) itself. That said, it is possible to show that the two equations above constitute a well defined initial-value problem. To do so, we differentiate the l.h.s. of Eq. (20) w.r.t. τ : The partial derivative ofq(τ, k 2 ⊥ ) can be obtained from the evolution equation in its integral form: =q(τ 0 , k 2 ⊥ ) +ᾱ s with τ s (k 2 ⊥ ) the inverse function of Q 2 s (τ ) that satisfies then τ s (Q 2 s (τ )) = τ . One can then compute the derivative w.r.t. k 2 ⊥ , evaluated at Q 2 s (τ ): Finally, we end up with the following differential equation for Q 2 s (τ ):  This equation is more convenient from a numerical point of view because one does not need to solve the implicit equation defining Q 2 s (τ ) at each small step in τ . Instead, one solves two coupled integro-differential equations of order 1.