Universality aspects of quantum corrections to transverse momentum broadening in QCD media

We study non-linear quantum corrections to transverse momentum broadening (TMB) of a fast parton propagating in dense QCD matter in the leading logarithmic approximation. These non-local corrections yield an anomalous super-diffusive behavior characterized by a heavy tailed distribution which is associated with L\'{e}vy random walks. Using a formal analogy with the physics of traveling waves, we show that at late times the transverse momentum distribution tends to a universal scaling regime. We derive analytic solutions in terms of an asymptotic expansion around the scaling limit for both fixed and running coupling. We note that our analytic approach yields a good agreement with the exact numerical solutions down to realistic values of medium length. Finally, we discuss the interplay between system size and energy dependence of the diffusion coefficient $\hat q$ and its connection with the gluon distribution function that is manifest at large transverse momentum transfer.

A Universal leading edge and front interior expansion with running coupling 50 1 Introduction The QCD jets that form in high energy scattering processes involving heavy nuclei such as heavy ion collisions (HIC) and electron-ion collisions undergo substantial final state effects caused by multiple interactions with hot or cold nuclear matter. An important experimental signature of the latter is the strong suppression of jets and high-p t hadrons, commonly referred to as "jet quenching" [1,2], that was observed in HIC at RHIC in the early 2000's [3,4] and more recently at the LHC [5][6][7]. In order to gain more insight on the transport properties of the quark gluon plasma (QGP) that is created in the aftermath of the collision, extensive studies of jet events and their substructure are currently carried out [8,9]. One of the key observable is the dijet acoplanarity that measures the amount of azimuthal decorrelation that a dijet system suffers because of final state interactions [10][11][12]. It results in particular in the broadening of the jet transverse momentum w.r.t. their initial direction of propagation. Although this effect remains illusive at the LHC due to the large contribution of multi-jet events that are formally resummed in the so-called Sudakov from factor [13], there are hints of momentum broadening at RHIC energies [13,14] (see also [15] for a similar discussion using substructure techniques). Transverse momentum broadening plays also an important indirect role in the process of radiative energy loss which is the main cause of the phenomenon of jet quenching.
Up until recently, transverse momentum broadening (TMB) was mostly studied at tree level. While multiple scatterings are resummed to all order in the limit where the in medium mean free path mfp is much smaller than the medium length L, the scattering rate is usually computed at leading order, i.e. α 2 s , in the limit of large momentum transfer q ⊥ T , where T is the plasma temperature. However, higher order corrections turned out to be enhanced by potentially large double logarithms of the medium length α s ln 2 (L) [16] and therefore must not be neglected. These corrections were shown to be related to a renormalization of the jet quenching parameterq ≡ d k 2 ⊥ typ /dt that measures the typical transverse momentum squared per unit time accumulate by the fast parton in the plasma [17,18]. On the other hand, they involve gluon fluctuations that may stretch beyond the in-medium correlation length all the way up to the medium length reflecting nonlocality of interactions beyond leading order. As a consequence, the diffusion coefficient is time dependent and results in an anomalous diffusion in transverse momentum space in contrast with normal diffusion at tree-level [19].
A systematic control of these logarithmically enhanced quantum corrections is not only crucial for precision phenomenology at RHIC and LHC, as well as at the EIC, but also for probing this novel quantum transport phenomenon in experiment. Some progress has been made in this direction. The single logarithmic corrections enhanced by α s ln (L) have been computed in the seminal paper by Liou, Mueller and Wu [16]. More recently, Zakharov pointed out that relaxing the soft gluon approximation, usually performed in the computation of the radiative corrections, may lead to non-negligible NLO corrections at RHIC or LHC kinematics [20]. Furthermore, the interplay between the dense and dilute limit has been studied in [21]. E. Iancu proposed a non-linear evolution equationà la JIMWLK to account for the dominant radiative corrections to all orders [18]. However, due to its complexity, there is, so far, no known solution (neither analytic nor numeric) to this evolution equation. Nevertheless, in the double logarithmic approximation (DLA), which aims at resumming the subset of these corrections whose general term takes the form α n s ln 2n (L), the equation simplifies and can be studied analytically or numerically both at fixed [16,17,22] and running coupling [23].
Another way to think about this resummation program is in terms of Wilsonian renormalization: the divergence of the quenching parameter (defined at a time scale τ 0 ∼ 1/T ) due to an additional soft and collinear gluon emission with lifetime between τ 0 and τ 0 + dτ is absorbed into a redefinition of the quenching parameter at a time scale τ = τ 0 + dτ [17]. This renormalization ofq is likely process-independent. For instance it has been shown that the leading, double logarithmic, radiative corrections to the medium-induced gluon spectrum of an off-shell quark can also be included by using the leading order BDMPS-Z spectrum [24][25][26] with the same effective (renormalized)q as the one involved in the radiative corrections to TMB [17,18]. Recently, the authors of [27,28] demonstrated that this universality also holds at single logarithmic accuracy.
In this work, we perform an analytic and numerical study of the large system size dependence of the jet quenching parameter and the medium saturation momentum that is related to the typical transverse momentum broadening. We derive novel analytic results that are grounded on the double logarithmic approximation for the quantum evolution ofq. The latter appears to suppress the sensitivity to the initial, tree-level value of the quenching parameter, a property that we shall refer to as "universality". We carried out the computation of all universal terms in the large L expansion of both the jet quenching parameter and the saturation scale.
The double logarithmic evolution equation is considered both in the fixed coupling approximation and with running coupling correction. In the former case, we detail the mathematical results discussed in our short letter [19]. More precisely, thanks to a deep connection between the fixed coupling DLA equation and the physics of the propagation of traveling waves into unstable states [29][30][31][32][33][34] -typically governed by the Fisher-Kolmogoroff-Petrovsky-Piscounoff (FKPP) equation [35,36] -, we were able to compute all the universal sub-asymptotic corrections borrowing analytic techniques developed in that context. While previous studies were limited to its linearized version, in our approach, we successfully obtained the asymptotic expansion for the full non-linear version of the evolution equation. As we shall show, the non-linearities yield a sizable corrections, especially when L is not asymptotically large.
These techniques are also successfully applied to the evolution equation for the quenching parameter with running coupling. We provide in particular a proof of the numerically motivated ansatz proposed by Iancu and Triantafyllopoulos in Ref. [23]. In that case too, we are able to push further the accuracy of the asymptotic expansion for large system sizes so that our formulas can be used even for realistic values of the in-medium jet path length.
We emphasize that the aforementioned system size dependence of the renormalizedq is different from the energy (or "x") dependence considered in theoretical studies ofq in the higher twist approach [37][38][39]. For high energy jets, with E qL 2 , the ln 2 (L) enhancement is the dominant contribution from radiative corrections, as we argue in Secs. 8.1-8.2.
The physical interpretation of our results has been largely discussed in [19]. We briefly summarize here our main findings, a more detailed discussion is provided in section 7. For asymptotically large system size, the transverse momentum broadening distribution reaches a universal distribution which satisfies geometric scaling. Namely, it is a function of k ⊥ /Q s only. This asymptotic distribution has interesting properties. First, its typical width that defines the saturation scale or equivalently the typical transverse momentum transfer (measured, for instance, by the median) scales like L 1/2+ √ᾱ s and thus, grows with the system size faster than L 1/2 . This is characteristic of a super-diffusive regime. Furthermore, the large k ⊥ tail of the asymptotic distribution exhibits a harder power spectrum 1/k 4−2 √ᾱ s ⊥ when compared to the Rutherford type, i.e. 1/k 4 ⊥ . These two propertiessuper-diffusion and heavy-tailed distribution-are reminiscent of Lévy distributions which describe the probability density for the position of a particle undergoing a Lévy flight. The self-similarity of gluon fluctuations and the non-local nature of the jet-medium interactions are two essential characteristics of Lévy flights, shared by other systems in various physical areas like optics [40,41], turbulence and polymer transport theory [42,43].
The article is organized as follows. In the first section we set up the formalism. We introduce in particular the notations and the evolution equations that we shall study in the rest of the paper. For readers who are only interested in the final results, the second section summarizes our main analytic findings regarding the quenching parameter and the saturation momentum as a function of the system size. In sections 4 and 5, we detail our analytic approach, based on the similarity with the problem of traveling waves, to compute the asymptotic and pre-asymptotic behavior of the quenching parameter and the saturation momentum for both the fixed and running coupling evolution. Section 6 presents an analytic approach to obtaining formulas approximating the exact solution down to realistic values of the system size. We briefly discuss qualitative aspects of the TMB distribution in section 7, and in particular the connection with statistical physics through the Lévy flight physical picture. Finally, in section 8, we comment on the relation between the operator definition of gluon distributions in DIS and quenching parameter and we discuss the energy dependence ofq at smaller jet energy. We highlight possible applications of our formula to phenomenology of heavy-ion collisions and small-x physics in our summary, Sec. 9.
cone direction 1 with a large longitudinal component of its momentum, P + ≡ E. More precisely, the momentum transferred from the medium is assumed to be small compared to the parton energy, that is, |k ⊥ | E. In this approximation, the transverse momentum broadening distribution of the high energy parton can be 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. [44][45][46][47]) by a Fourier transform as follows and is a path ordered Wilson line along x + and t a are the SU (3) generators in the fundamental representation if the fast parton is a quark. A Feynman graph representation of this treelevel calculation is shown in Fig. 1a. The medium background field A − has support in [0, L] where L is the length of the medium and ... stands for an ensemble average over the medium color charge configurations. Note that Eq. (2.1) may be encountered in a variety of high energy scattering processes. For instance, in gluon saturation physics at small-x probed in inclusive DIS or forward particle production in high energy collisions such in pA collisions for instance, A − is the classical field produced by a boosted nuclear/proton target. It also arises in jet production in heavy ion collisions. In that case, A − is the classical field generated by the color charges of the QGP.

Leading order: independent multiple scattering approximation
In the independent multiple scattering approximation the background field correlations are Gaussian [46]: Here, the leading order collision rate C encodes the Coulomb tail where n stands for the density of scattering centers. The expression (2.5) for the collision rate is valid in the perturbative regime for q ⊥ µ. Throughout this paper, µ denotes a generic non-perturbative scale. For a weakly coupled quark-gluon plasma in thermal 1 We use light cone coordinates defined as Figure 1: (Left) An illustration of a tree-level diagram [48] contributing to Transverse momentum broadening (TMB), where vertical gluons depict multiple scatterings between the incoming quark and the medium scattering centers. The dotted vertical line represents the cut between the amplitude and its complex conjugate (c.c.). (Right) Typical NLO correction to the tree-level distribution due to the radiation of a real gluon of light-cone energy ω and transverse momentum k ⊥ . The transverse vector x ⊥ is the difference in transverse position of the quark in the amplitude and its c.c. It is the Fourier conjugate of the measured transverse momentum broadening.
equilibrium at temperature T , µ is of order the Debye mass m D ∼ gT . On the other hand, in the case of a boosted nucleus, µ is a scale of order the inverse nucleon size. The precise form of the collision rate in the non-perturbative domain will not be discussed in the present article. We will mainly focus on the TMB distribution in the perturbative regime where the effect of leading quantum corrections is to strongly suppress non-perturbative physics, resulting in a universal behavior as shall be discussed in detail. Under the independent multiple scattering approximation, S(x ⊥ ) exponentiates as where the leading order diffusion coefficientq (0) readŝ up to powers of µ/Q suppressed terms. The color representation R = A, F of the hard parton enters through the Casimir C R . Thereafter, we shall write the tree-level quenching parameter in a more generic way asq (0) =q 0 ln(Q 2 /µ 2 ), with the constantq 0 and µ being model dependent. For instance, for a weakly coupled QGP the bare quenching parameter q 0 and the infrared transverse scale µ 2 can be obtained from the hard thermal loop value of the collision rate [49]. They read respectivelyq 0 = α s C R m 2 D T and µ = m D e −1+γ E /2 [50], with T the plasma temperature, m D the Debye mass.
Eq. (2.6) encompasses two regimes separated by an emergent scale, the saturation scale Q 2 s (L) (in analogy with gluon saturation at small-x): the regime k ⊥ Q s is characterized by rare hard events that correspond to S 1. The TMB distribution falls like a power law, P(k ⊥ ) ∝ 1/k 4 ⊥ that is characteristic of Rutherford-like scattering. When k ⊥ Q s the unitarity bound is saturated, i.e., S 1 and all scattering centers contribute equally. For k ⊥ ∼ Q s , the k ⊥ distribution is typically Gaussian, P(k ⊥ ) ∝ exp(−k 2 ⊥ /Q 2 s ). As is customary in small-x physics [51,52] or in Molière theory of multiple scattering [50,53,54], the saturation scale is mathematically defined by the relation S(x 2 ⊥ = 1/Q 2 s (L)) ≡ e −1/4 , or equivalently, Eq. (2.7) can be solved for moderate values of L iteratively and one finds approximately Q 2 s ∼q 0 L ln(q 0 L/µ 2 ) at tree level for the saturation scale. Notice that the saturation scale Q s depends implicitly on the color representation of the high energy parton via the Casimir dependence of the quenching parameter.

One-loop corrections and the double logarithmic resummation
TMB at one loop. Let us now discuss the transverse momentum distribution at one loop order. At order α s , one must include two contributions, one real and one virtual gluon attached to the incoming parton. A typical Feynman diagram representation of a real NLO correction that contributes toq is presented in Fig. 1b. The limit in which the additional gluon is soft was addressed in detail in [17] (see also the appendix A of [16] for a calculation in Zakharov's formalims). Leaving technical details aside for the sake of clarity, we may write schematically: where P (0) (k ⊥ ) is the tree-level TMB distribution given by Eqs. (2.1)-(2.6). Even though Eq. (2.9) looks like a standard perturbative expansion, one should keep in mind that each term actually resums to all orders powers of α 2 s nL due to multiple scatterings. One can show that the leading contribution from α s P (1) (k ⊥ ) comes from radiative corrections which are quasi-local, in the sense that they can occur everywhere inside the medium over the path length L of the high energy parton. The average k 2 ⊥ associated with such contributions reads where ω and τ are respectively the energy and the lifetime of the gluon fluctuation. This lifetime is bounded from above by the typical formation time of a medium-induced emission triggered by multiple soft scatterings. The latter bound ensures that the radiative process is triggered by a single scattering with medium constituents. The other boundaries of the double integral will be specified below. This double integral has a typical double logarithmic structure, k 2 ⊥ ∼ α s ln 2 , and therefore, the smallness of α s can be compensated by the large logarithm squared, spoiling the convergence of the series (2.9). Note also that k 2 ⊥ is enhanced by the system size L, as a result of the quasi-locality of the double logarithmic quantum corrections [17].
In this paper, we are mainly interested in the resummation of double logarithmic contributions of the form (2.10) to all orders in the series (2.9). This resummation can be performed at the level of the quenching parameter itself, or in other words, the double logarithmic corrections exponentiate owing to the fact that the logarithmic τ integral is dominated by the regime τ L (for more details see discussion in [17]): withq (1) ∼q 0 ln 2 ,q (n) ∼q 0 ln 2n , etc. A complete NLO computation of the TMB distribution would require both the resummation of the double and single logarithms via this exponentiation property, and a proper matching with the fixed order result (2.9). Such calculation is beyond the scope of this paper. We refer the interested reader to Refs. [16,27] where some aspects of the single logarithmic contribution to momentum broadening and radiative energy loss are discussed.
Double logarithmic phase space forq. We now detail the double logarithmic phase space for the quenching parameterq. In order to specify these boundaries, it is more convenient to perform the change of variable ω → k 2 ⊥ = 2ω/τ . From the logarithmic corrections, the quenching parameter acquires a τ and k 2 ⊥ = 1/x 2 ⊥ dependence. At one loop, the double logarithmic contribution readŝ Here τ 0 is a cut-off time scale of the order of the mean free path, which reflects the uncertainty of this calculation due to the non-perturbative physics of the plasma. Q s (τ ) is the saturation scale at the lifetime τ , defined through the relation (2.8) or in terms of the functionq(τ, k 2 ⊥ ): It is instructive to estimate this integral using a constant tree-levelq value, namelyq(τ, k 2 ⊥ ) = q 0 leading to Q 2 s (τ ) q 0 τ . The consequences of this approximation on the evolution of the quenching parameter will be discuss in details in this paper. The single hard scattering condition τ ω/q becomes k 2 ⊥ Q 2 s (τ ) with our new variables. Note that the color factor is N c since we are dealing with the scattering of a gluon at the one loop order.
The main difference with the standard double log encountered in DGLAP is that the collinear and soft logs talk to each other through the saturation line that plays the role of physical cutoff for the collinear singularity as a result of coherence effects from multiple scattering.
If one is interested in k ⊥ Q s (L) we would have (neglecting the logarithmic dependence of theq (0) ≡q 0 at leading order for simplicity) with µ 2 =q 0 τ 0 andᾱ s = α s N c /π, where the microscopic scale τ 0 is related to the in-medium mean-free-path which for a thermal plasma scales as (g 2 T ) −1 at weak coupling. This double integration corresponds to the area of the right trapezoid depicted in Figure 3 (left panel). When the upper limit of the k ⊥ integration falls below Q 2 s , that is, k 2 ⊥ < Q s (L) 2 , one is left with the area of a triangle and one obtains the following double log (cf. Figure 3 (right panel)) To obtain the corrections to the typical value of transverse momentum broadening we must evaluateq at k 2 ⊥ = Q 2 s (L) q (0) L which yields the Liou-Mueller-Wu result [16] k 2 Resummation. In the double logarithmic accuracy (DLA) these radiative corrections can be resummed to all orders via an evolution equation ordered in τ [16][17][18]: whereq (0) (τ 0 , k ⊥ ) corresponds to the tree-level initial condition. The strong coupling constant appears inside the integral over k ⊥ to account for its running with the transverse scale. It is convenient to re-express these two equations in terms of the logarithmic variables Y = ln τ τ 0 and ρ = ln This non-linear evolution equation resums the double logarithms α s Y ρ to all orders. Also, it is valid in the large N c limit which is reflected in the overall N c factor absorbed in the constantᾱ s . In principle, since the definition of Q s is "flavor" dependent, there should be a coupling between the evolution of the quenching parameterq F in the fundamental representation and the adjoint oneq A . In this paper, we do not consider the effect of such a coupling (which is beyond DLA) and focus on the evolution ofq A only. At this accuracy, the fundamentalq can be obtained fromq A usingq F = C F /C AqA . A graphical representation of this evolution equation is displayed in Figure 2.
For the running coupling evolution, we use the one loop beta function to determine the ρ dependence ofᾱ s :ᾱ ... Figure 2: An illustration of multiple radiative corrections considered in this paper. Each block represents a tower of gluon fluctuations triggered by a single scattering with a medium constituent, with strongly decreasing lifetime and transverse momentum along the cascade (so the transverse size of a gluon increases from the parent to its daughter). The exponentiation resums several such blocks over the path length L of the incoming hard effective dipole with transverse size 23) where N f = 5 is the number of quark flavors.
Dense and dilute regime. In determining the transverse momentum distribution we have to distinguish between the dense regime ρ < ρ s (Y ) and the dilute one ρ ρ s (Y ). In the latter, recall that the quenching parameter is function of two independent variables Y = ln(L/τ 0 ) and ρ. In the dense regime, however, there is a subtlety in the choice of variables. Given a general solutionq(Y, ρ), the variable Y is no longer an independent function of ρ. Again, this variable is related to the upper limit of the τ integral in the DL phase-space as illustrated in Figure 3. Therefore, the logarithmic variable Y must be fixed such that the quenching parameterq which appears inside the forward scattering amplitude is only a function of ρ. The relevant time scale τ s (k 2 ⊥ ), or in logarithmic variables Y s (ρ), at which the quenching parameterq must be evaluated is the largest time allowed by the saturation condition: The function Y s (ρ) is then the inverse function of ρ s (Y ). Again, in the dilute regime, the value of Y is fixed by the typical path length L of the hard parton inside the dense medium [21], This is illustrated in Figure 3 where we see that for k 2 is defined by τ s = k 2 ⊥ /q(τ s ), whereas for k 2 ⊥ < Q 2 s (L) we simply have τ < L.
In summary, the function of k 2 ⊥ (or ρ) to be used in the forward dipole amplitude is: As we shall see in Sec. 6, this function is continuous and derivable in ρ = ρ s (Y ), but the second derivative is not in general continuous.
Linearization. Exact analytic solutions of the non-linear system (2.17)-(2.18) are in general difficult to obtain. However, there are analytic solutions to the fixed coupling, linearized problem that consists in approximating Q 2 s (τ ) q 0 τ in the lower bound of the k ⊥ integral in Eq. (2.17) [23,55]. Under these approximations, the integro-differential equation forq decouples from the implicit equation satisfied by Q s . The saturation scale is defined instead by In terms of the logarithmic variables Y and ρ,q(Y, ρ) satisfieŝ The analytic solutions forq(Y, ρ) can be obtained by iterations [16,23,55]. For constant tree-level initial conditions,q (0) (ρ) =q 0 , the solution readŝ while for leading twist tree-level initial conditionsq (0) (ρ) =q 0 ρ, one obtainŝ where the I n (x)'s are the modified Bessel functions of the first kind. Note that since we consider the linearized evolution equation, the solution with initial conditionq (0) (ρ) = q 0 ρ+const. is simply the linear superposition of the two solutions above. In the asymptotic limit, whenᾱ s Y ρ 1 the quenching parameter behaves roughly likê A detailed discussion of the asymptotics of these solutions will be presented in Section 4.1. Figure 3: The double logarithmic phase space for gluons fluctuations with lifetime τ and transverse momentum k ⊥ in the dilute (left) and dense (right) regimes. The choice between these two regimes depends on the final momentum k ⊥ of the double logarithmic cascade. For the TMB distribution, we set k ⊥ = 1/x ⊥ where x ⊥ is the initial effective dipole size.
In this sketch, we have simplified the saturation line Q s (τ ) by its tree-level formq 0 τ .
3 Asymptotics of TMB: highlights of the main results for E qL 2 This section is meant for the reader who is less interested in the formal derivations detailed in the next two sections, but rather in the analytic formulas we have obtained that can be used for heavy-ion or small-x phenomenology. We will focus on high energy limit: E ω c ≡qL 2 and postpone the discussion of the opposite regime in Section 8.2.
In the following we summarize our analytic results for the asymptotic expansion of the quenching parameter and the saturation scale in both fixed coupling and running coupling scenarios. We display our results forq(L, k 2 ⊥ = 1/x 2 ⊥ ) as a function of ρ = ln(k 2 ⊥ /(q 0 τ 0 )) and Y = ln(L/τ 0 ). Similarly, the asymptotic expansion of the saturation scale is written in terms of ρ s (Y ) = ln(Q 2 s (L)/µ) so that

Fixed coupling
For the fixed coupling evolution, the saturation scale is given by where c = 1+2 ᾱ s +ᾱ 2 s +2ᾱ s is the celerity of the front and κ is a non-universal integration constant to be determined numerically. Defining

Running coupling
In the running coupling case, we present analytic results that are relevant for realistic values of Y (not too large, typically smaller than Y = 5). For such values, the effects of the non-linearity are mild and one can focus on the linear evolution. The saturation momentum reads in that case: The quenching parametersq(Y, x) which enters inside the TMB distribution is given bŷ Note that compared to the fixed coupling result, we have not expanded in powers of Y each coefficient of the polynomial functions in x which appear inq(Y, x). The reason is that the convergence of the development of ρ s is slower (due to the 1/6 power instead of 1/2). Therefore, one must keep the entire ρ s development (3.5) in order to reach realistic values of Y . With this expression, we are able to accurately estimate the k ⊥ -distribution over 4 orders of magnitude around the saturation momentum, as shown Figure.

Geometric scaling at fixed coupling
In the regime where the energy E of the propagating parton is much larger than the characteristic frequency ω c =qL 2 , the broadening distribution does not depend on E. As a result, the physics will be governed by a single scale Q s (L) that is a function of system size and local medium properties such as the temperature T . Furthermore, the k ⊥ -distribution will depend asymptotically on a single scaling variable k ⊥ /Q s over a wide range of k ⊥ 's. The latter property will be refer to as geometric scaling.
To study the asymptotics of this regime, we shall take the formal limit E → ∞ to compute the asymptotic behavior of the quenching parameterq(L, k 2 ⊥ ) and the saturation scale Q s (L) for large system sizes. Our calculation will be accurate up to terms of order 1/ ln(L). Solving exactly the non-linear evolutions forq is a difficult task. Therefore, in order to gain insight on the generic properties of solutions we shall first analyze the linearized version of the evolution equation for which exact analytic results can be obtain for arbitrary system sizes. This exercise will serve as an anchor point for the non-linear case where only asymptotic results will be derived.
Furthermore, we will take advantage of a mathematical equivalence between the evolution of the quenching parameter and the formation of traveling wave-front in non-linear physics to discuss the onset of the extended geometric scaling.
In section we focus on the fixed coupling case. The running coupling will lead to a different scaling which will be discuss in detail in Section 5.

Asymptotics of the analytic solution to the linearized equation
We remind the reader that the linearized evolution equation is obtained by replacing Q s in the lower bound of the k ⊥ integral in r.h.s. of Eq. (2.17) by its leading order valueq 0 τ . This turns out to be a good approximation since non-linear effects do not alter the leading asymptotic behavior as we shall see.
Thanks to the analytic solutions (2.29) and (2.30), we want to show that the linearized evolution equation satisfies the following scaling property at large Y : for some function ρ s (Y ) to be determined. Using the asymptotic expansion of the Bessel functions, In order to get the correct scaling function when Y → ∞, it is necessary to keep track of the sub-leading terms in the expansion (4.2). Replacing ρ = x + Y and expanding for large Note that one obtains the same result starting from Eq. (2.30), the difference between the two initial conditions is encoded in the constant term in ρ s (Y ). This is an important observation as it illustrates the universality of the asymptotic solution forq(Y, ρ), namely, the fact that it loses sensitivity to the tree-level initial condition at large Y = ln L/τ 0 . One can even infer from the exact analytic solutions for Q 2 s (Y ) the order of the non-universal coefficients. For a constant initial condition, we have (using Eq. (2.27)): while for the initial conditionq (0) (ρ) =q 0 ρ, we get We conclude then that for the linearized evolution equation, only the first two terms are universal, whereas the constant and O(Y −1 ) terms depend on the initial condition. As we shall see, these statements remain correct in the non-linear case, with two important differences though: the coefficients of the Y and ln(Y ) terms are modified and another universal power 1/ √ Y appears in the development. Finally, we emphasize that not only does Eq. (4.4) account for the scaling limit, given by but it also encompasses the sub-asymptotic corrections via the 1/Y suppressed term inside the exponential and the prefactor. The latter informs us about the range of geometric scaling, that is, x Y ∼ ρ s , or in terms of physical variables k 2 ⊥ Q 4 s /µ 2 . Or course, in this analysis the scaling is not exact since the scaling variable is x = ρ−Y which involves Y rather than ρ s . We shall see that in the non-linear equation Y is replaced by ρ s (Y ) resulting in a scaling with ρ s , i.e. Q s , only .

Asymptotic behaviour from non-linear wavefront formation
In this section, we extend the previous results to the non-linear evolution equation for q(Y, ρ), exploiting a formal analogy with the physics of traveling waves propagation into unstable states [29][30][31][32][33][34]. Such a mathematical connection turned out to be also fruitful for the study of the asymptotic behaviour of the saturation scale resulting from the non-linear BK evolution [56][57][58][59][60].
The study of the linearized evolution equation suggests an asymptotic scaling solution of the formq where f is only a function of the scaling variable The connection to traveling wave physics is made more transparent with the dipole S-matrix defined in Eq. (2.2): (4.14) where ρ = − ln(x 2 ⊥q 0 τ 0 ). For ρ ρ s , S saturates at 1 owing to the unitarity constraint and drops to zero when ρ ρ s as illustrated in Figure 5. Hence, as Y increases ρ(Y )  Figure 5: Illustration of the traveling wave and its front propagation. The front propagates from the left to the right as Y increases. The domains studied analytically in this paper are represented: the interior of the front, corresponding to ρ − ρ s (Y ) Y α , where the non-linearities become important, and the leading edge regime with ρ − ρ(Y ) ∼ Y α . In the fixed coupling evolution, α = 1/2 while for the running coupling one, α = 1/6. increases and the evolution of S as function of Y can be interpreted as the propagation of a front along the ρ axis with speedρ s (Y ) = const. and Y playing the role of the time. In the scaling regime, the wave propagates to the right while retaining its (universal) shape. The calculation of the deviation with respect to this uniformly translating profile at asymptotic times Y can be done using techniques borrowed from non-linear physics of wavefront formation.

Traveling wave solutions
We return now to the full non-linear problem for which Eq. (2.28) is replaced bŷ The difference with Eq. (2.28) lies in the lower bound for the ρ integral, where we substituted Y → ρ s (Y ). In order to proof the existence of a family of front solution we insert the scaling ansatz (4.11) into the fixed coupling evolution equation Eq. (4.15) forq(Y, ρ). We obtain the following integro-differential equation which upon differentiation w.r.t. x reduces to a second order differential equation: Now, requiring thatq satisfies the scaling relation (4.11) implies that the Y dependent coefficients must be constant:ρ for large enough Y . Thus, we need to solve a the second order differential equation: The basis of solutions are given by exponential functions of the form exp(βx) with β such that . (4.20) Rejecting imaginary values for β (non-oscillatory solutions) provides the constraint c 1+2 ᾱ s +ᾱ 2 s +2ᾱ s on the front velocity, which is equivalent to β β c = −ᾱ s + ᾱ s +ᾱ 2 s . The value of c dynamically selected by the system depends on the initial condition at Y = 0. It is straightforward to see that for asymptotic solutions of the form exp(βx), the function 1 − S(0, ρ) behaves like exp((β − 1)ρ) at large ρ (note that β − 1 < 0 forᾱ s > 0). If the initial condition forq is such that 1 − S(0, ρ) decays faster than exp((β c − 1)ρ), i.e. q(0, ρ)e −βcρ → 0 for large ρ, then the value of c chosen by the system is the minimal one [30,31,33,34,[61][62][63]68]: From now on we shall drop the subscript c of β c since we will only consider physical initial conditions that satisfy this requirement. For this critical value c, the solution f is given by f (x) = e βx (a 1 + a 2 x) with a 1 and a 2 two integration constants. They are fixed by the very definition of ρ s , which is given by the relation that implies that f (0) = 1 and f (0) = (c−1)/c on the saturation line x = 0, i.e. ρ = ρ s (Y ) yielding These results should be contrasted with those obtained for the linear evolution equation and reported in Eq. (4.10). First, the non-linear evolution imposes a single scale in the scaling limit, Q s . Therefore, the function f is a function of In terms of physical variables, the scaling variable is k 2 ⊥ /Q 2 s in the non-linear case and k 2 ⊥ /(q 0 L) in the linear case. Interestingly though, the shape of the scaling function f and the sub-asymptotic corrections toq(Y, ρ) are the same within the double-logarithmic approximation in which c 1 + 2 √ᾱ s and β √ᾱ s given that we work in the weak coupling limit, namely,ᾱ s 1 .

Logarithmic shift of the front: study of the leading edge domain
The traveling wave solution f describes a uniformly translating front invading the large ρ domain where 1 − S vanishes. The corrections to the velocityρ s for non-asymptotic values of Y are driven by the leading edge domain of the solution, that is the region where the front forms as a result of the growth of perturbations around the "unstable state" 1−S = 0.
To obtain the behavior of the front on the leading edge, we look for a solution of the form [57,64]q This ansatz will allow us to describe the leading edge where z = x/Y α ∼ 1. Assuming α to be positive and since Y is large, it follows that x = ρ − ρ s is large as well. Thus, we are indeed focusing on the front formation region. The power α describes a diffusion-like spreading of the traveling wave in the leading edge domain. Plugging this ansatz inside the differential equation satisfied byq and using the definitions of c and β given by Eq. (4.21) in terms ofᾱ s , one gets the following differential equation for G α : The homogeneity constraint of this equation implies that α−1 = −α and thus α = 1/2 anḋ The general solution of this equation can be expressed in terms of confluent hypergeometric functions 1 F 1 (a; b; x) and Hermite polynomials H(n; x): (4.28) The unknown constants δ 1 , a 1 and a 2 can be determined from the boundary conditions at z = 0 and z → ∞. G has to decay faster than a power law at large z [64], meaning that the constant a 2 must vanish as it multiplies the exponentially growing solution 1 F 1 (a; 1/2; x) ∼ e x x a−1/2 . The initial condition at z = 0 is obtained by requiring that the solution G matches with the scaling limit f for large x, so that The small z behavior of the Hermite function is In order for the constant term to vanish, the argument of the Γ function must be a negative integer, constraining δ 1 to be of the form δ 1 = − 3+2n 1−β with n a non-negative integer. Finally, since G must always be positive (no nodes) [64,68], this fixes n = 0 and therefore This yields the following sub-asymptotic correction for the shape and position of the wavefront: As for the case of the pulled front propagation into unstable states, our analysis of the leading edge domain shows that for non-asymptotic "times" Y , the position of the wave-front undergoes a logarithmic shift, whose coefficient is given by the constant δ 1 . We observed the same logarithmic sub-asymptotic corrections in the linearized evolution. However, in the latter the coefficient was −3/2, whereas it is of order −3/2(1 + √ᾱ s ) in the non-linear case. As argued in [19], this is parametrically larger than single logarithmic corrections that we neglect in this paper. On the other hand, comparing Eq.

Leading edge vs. front interior expansion
So far, we have obtained two distinct asymptotic solutions to the non-linear evolution equation forq(Y, ρ): f and G. The former is valid in the interior of the wave front corresponding to the saturation regime where S ∼ 1 and 0 < x 1, whereas, the latter applies to the leading edge domain where x ∼ Y 1/2 1 and S 1. This solution reads, in terms of x and Y (that is, in the front rest frame): Expanding for x 2 Y (or equivalently z 1), one findŝ This suggests the existence of an other expansion scheme of the form . This expansion is dubbed the "front interior expansion" as it focuses on the behaviour of the solution in the rest frame of the front and near the saturation region (z 1) [64]. In this region, it is not allowed to linearize the evolution equation.
On the other hand, one can also study the corrections to the solution G in the leading edge domain (z ∼ 1). The natural expansion scheme there will be referred to as "leading edge expansion" [64], and it readŝ with G −1 ≡ G. Each term in this series resums to all orders powers of x/ √ Y . By expanding in power series these analytic functions of z, one notices a relation between the leading edge and the front interior expansion. This is represented in Table 1: each function in the leading edge series accounts for a diagonal in the infinite triangular matrix of coefficients f j n . The function G −1 resums the first red diagonal, the function G 0 resums the second diagonal and so forth. There is a systematic way to compute each term in the leading edge, front interior andρ s expansion. As we shall see, these two expansions are related by matching conditions and the leading edge expansion constrains the coefficients of the asymptotic expansion of the saturation scale ρ s .
We already know f 0 and G −1 . The function f 1 and f 2 can be obtained without difficulty by plugging the front interior expansion in the evolution equation forq and usingρ s c + δ 1 /Y . One gets the following differential equations: The initial conditions for f 1 and f 2 are given by the definition of Q s which leads Some features of this calculation are generic to all f n functions. First, they all satisfy a second order differential equation of the form f n (x) = ... whose r.h.s. is determined by the f i (x) with i < n. The initial conditions are always provided by the definition of the saturation scale which yields f n (0) = 0 for n ≥ 1. Therefore, all the f n (x) are polynomial functions of x, with degree at most n + 1.
To better understand the interplay between the leading edge and front-interior expansion, we compute the functions G 0 and the next term in the development ofρ s . Inserting the leading edge expansion into the evolution equation forq and the development ofρ s = c+σ s , one finds a differential equation for G 0 whose homogeneity condition constrains the form ofσ s to beσ The homogeneous part of the differential equation satisfied by G 0 is similar to the one satisfied by G −1 , while its inhomogeneous term depends on G −1 and its derivatives: To solve this differential equation, it is convenient to perform the change of variable Replacing the known values of β, δ 1 and G −1 in terms of c, one finds the following differential equation for g 0 (u): and C = √ c − 1. The integration constants are fixed by matching the leading edge expansion with the front interior. More concretely, since f 0 (x) = 1 + O(x) and f 1 (x) = 0, the development at small u of g 0 must be g 0 (u) = 1 + O(u) (terms of order u 1/2 are prohibited). The constant δ 2 is fixed by demanding that at large u, g 0 (u) diverges no slower than e u u −3/2 [64]. We then obtain for δ 2 the value while the function g 0 can be expressed in terms of hypergeometric functions 2 F 2 and 1 F 1 [58,64]: It is enlightening to expand the function G 0 (z) in powers of z. From the expression above, one finds Combined with Eq. (4.35), one observes that the function G 1/2 and G 0 resums respectively the leading and sub-leading powers of the polynomial functions f n (x) for all values of n, as illustrated in Table 1. In Fig. 6, our analytic results are compared to numerical simulations of the non-linear system (4.15). On the left figure 6a, the dashed blue line corresponds to the dipole Smatrix in the rest frame of the front, i.e. as a function of k 2 ⊥ /Q 2 s (L), given by the analytic expression ofq including the first two terms in the leading edge development, while the red P P P P P P P P P f n (x)  Table 1: Relation between the front interior and leading edge expansion. The function G −1 resums the red powers in x to all orders, which correspond to leading powers of the f n polynomials, the function G 0 resums the blue terms, the next-to-leading powers, and so forth. curve is the numerical solution. Even for this rather realistic value of Y = 5, the leading edge development shows a rapid convergence (the analytic curve could be further improved around the transition at x ≈ 0 by including more terms in the front interior expansion, see the discussion in Sec. 6). The black dotted line is the asymptotic scaling limit given by the function f . The right plot 6b displays the velocity of the frontρ s as a function of time Y , as given by our asymptotic expansioṅ and its truncation up to order Y and 1/Y , compared to the numerical solution. Once again, the convergence of the development is very good, down to small values of Y ≈ 2 ÷ 5.

Modified geometric scaling with running coupling
In this section, we consider the double logarithmic evolution of the quenching parameter including the running of the strong coupling constant. One would naively expect that running coupling corrections modify the fixed coupling asymptotic results for ρ s (Y ) via terms of order α s Y (single logarithmic corrections instead of double logarithmic ones). However, it is not the case as the power structure of the asymptotic development of ρ s at large Y is dramatically modified compared to the fixed coupling scenario. This is mainly a consequence of the evolution equation with running coupling belonging to a different universality class than the FKPP equation of the fixed coupling evolution.

Proof of the Iancu-Triantafyllopoulos's conjecture
In this subsection, we derive the scaling limit and the sub-asymptotic corrections of the solution to the linearized evolution equation with running coupling and thus, provide the proof of the expansion for ρ s (Y ) conjectured in [23] based on a numerical analysis. The equation we aim at solving is as follows where we have approximated the running coupling bȳ Although neglecting ρ 0 is justified asymptotically by the fact that for large Y , ρ > ρ s ρ 0 the following analysis can be extended easily by simply shifting ρ → ρ+ρ 0 and Y → Y +ρ 0 . Following [23], we introduce a new variable u such that u = ln(ρ/Y ) and f (Y, u) =q(Y, ρ) .
We will see that the g n functions are polynomial functions of z of degree n. After the Plugging Eq. (5.4) in Eq. (5.5) and differentiating with respect to u, we find Since g 0 (u) = 1 is a polynomial, it is straightforward to show by recursion that the g n (u) functions are polynomial in u, assuming that they do not diverge exponentially at large u. They are also positive for u ≥ 0. On the tree-level saturation line, that is, u = 0, the differential equation satisfied by g n implies This equation is very interesting: it shows that the n + 1 coefficient of the saturation line is related to a Laplace transform of the polynomial that multiplies the previous power Y n . Hence, if we can obtain the asymptotic form of g n (u) when n → ∞ we then may use the steepest descent method to integrate over u. This is the method that we shall follow to derive the scaling limit ofq(Y, ρ).

Scaling solution
First, we would like to have a visual insight on the shape of the polynomial functions g n (u). The functions g n (u) for 0 ≤ n ≤ 5 are shown in Fig. 7a. They do not seem to follow any organizing principle. However, once one considers the rescaling g n (u) → e −u g n (u/n)/g n (0), inspired by Eq. (5.7), we observe, as shown in Fig. 7b, that the corresponding curves all tend to lie on the same universal scaling limit as n becomes large. It is then tempting to conjecture that there exists a scaling behavior of these polynomial functions at large n. Thus, inserting the following ansatz: with the new variable ν ≡ nu.
-25 -Assuming a separation of variables we can write: g n+1 (0) = c (n + 1) 2 g n (0) , (5.10) and where c is a constant de be determined. Comparing with [23], we should expect c = 4. The solution that satisfies h(0) = h (0) = 1 reads In fact, c = 4 appears to be a critical value (it is a double root) that separates oscillatory solutions from exponentially decaying ones. Like with fixed coupling evolution, for realistic initial conditions, the value of c chosen by the system at large time is the smallest among all physical values (those which do no lead to an oscillatory behaviour). We shall see in the next subsection that the traveling wave interpretation of the running coupling evolution allows us to understand this criterion in the same way as at fixed coupling. Hence, at large n we readily find g n (0) = 4 n (n!) 2 . (5.13) We are now equipped to derive the leading asymptotic behavior of the quenching parameter from the expansion (5.4). Using Eq. (5.18) Expressing the above result in terms of the variable x = ρ − Y by approximating z = ln(ρ/Y ) x/Y 1, we find We recognize in the first factor the saturation scalê Turning now to the exact result forq(Y, ρ), one obtainŝ where the second term inside the parenthesis comes from the nu/2 term in Eq. (5.14).
Not surprisingly, we recover the fixed coupling exponent β ∼ √ᾱ s ∼ b 0 /Y . However, contrary to the fixed coupling case, the scaling variable is no longer x, but x/ √ Y . The leading behavior of the saturation scale, which appears as a pre-factor in (5.21) is

Sub-asymptotic corrections
In order to draw a complete picture of the solution we need to address the sub-asymptotic terms. Guided by the scaling analysis of the previous section, we may look for a solution of the form The pre-factor in the r.h.s. encodes the pre-asymptotic correction to g n (0). Now, we define the functionȧ(n) as g n (0)n 2 4g n−1 (0) ≡ eȧ . (5.25) Using the dotted notation for the discrete derivative with respect to n, we havė a(n) =L n (0) + ln n 2 − ln 4 . with L n (0) = ln g n (0). In the scaling limit, we haveȧ = 0 by construction, thereforeȧ quantifies the deviation of the ratio ( Similarly to the fixed coupling pre-asymptotic analysis, we may assume a diffusion-like Ansatz h(ν, n) = n α G y = ν n α . (5.28) It follows thatḣ where the second term can be neglected at large n since we expect α > 0. A similar argument allows us to neglect other G and G terms in what follows. Hence, G obeys the equation The homogeneity condition implies that α = 1/3, eȧ 1 +ȧ = 1 + βn −2/3 , and G (y) = 1 4 The exponent of the diffusive ansatz (5.28) is interesting as it differs from the fixed coupling evolution equation. Indeed, in terms of the scaling variable ν and Y , the function G is a function of ν/n α ∼ ν/Y 1/6 . This 1/6 power, smaller the power 1/2 at fixed coupling, changes the structure of the asymptotic expansion of ρ s quite dramatically. Turning back to the equation (5.32), the latter is the Airy equation with a regular solution at large y that reads In addition, in order to match onto the scaling solution we must impose that G(y) ∼ y at large n. This condition fixes the value of β to be where ξ 1 = −2.338... is the rightmost zero of the Airy function.
In terms of the original variables we finally obtain at large n g n (z) ∝ g n (0) n 1/3 Ai ξ 1 + 1 2 n 2/3 u , (5.35) with d dn ln g n (0) = −2 ln n + ln 4 + |ξ 1 |n −2/3 , (5. 36) or upon integration ln g n (0) = −2n(ln n − 1) + n ln 4 + 3|ξ 1 |n 1/3 + const. , (5.37) which is the form conjectured by the authors of [23] to fit their numerical results. Using again the steepest descent method as in the previous subsection, one obtains the subasymptotic corrections to ρ s (Y ): The analysis of the large n behaviour of the polynomial functions g n (u) enables us to obtain the first sub-asymptotic correction to the saturation scale, in agreement with the numerical findings in [23]. One could follow the same method to obtain the next corrections in the development of ρ s . Nevertheless, the traveling wave method detailed in following subsection turns out to be more efficient, since it offers a systematic way to compute the sub-asymptotic corrections and it can be extended to the non-linear case in which the saturation boundary in the evolution equation is fixed at ρ s (Y ) instead of Y .

Asymptotic analysis of the non-linear equation
We now consider the non-linear evolution equation with running coupling Even though the running coupling evolution of the quenching parameter does not fall into the same universality class as the fixed coupling equation (or the FKPP equation), we will follow the same strategy as in section 4.2 in order to get the scaling limit and the sub-asymptotic corrections.

Scaling limit and diffusive deviations
Modified scaling variable. We start by re-deriving the scaling limit given by Eq. (5.21) for the linear evolution equation. Inspired by the form of the scaling variable (ρ − Y )/ √ Y in that case, we try the following scaling form This scaling variable χ should be contrasted with the corresponding variable x = ρ − ρ s (Y ) in the fixed coupling evolution. Inserting this ansatz into the non-linear evolution equation forq(Y, ρ) with running coupling, one finds that the function f satisfies Setting χ = 0, and using the definition of the saturation momentum f (0) = 1, the velocity of the front behaves asρ Differentiating once more w.r.t. χ and using the asymptotic behaviour ofρ s given above yields with c = f (0). Expanding in powers of Y on both side, and extracting the leading 1/Y 1/2 terms, one gets As in the fixed coupling case, the value of c is fixed by requiring the discriminant of this differential equation to vanish, so that We have then recovered the results (5.21), namely, Diffusive ansatz around the scaling limit. Following the method detailed in Sec. 4.2, we need to compute the perturbations in the leading edge domain of the front in order to get the next term in the development of ρ s at large Y . We consider the following diffusive ansatz, similar to the fixed coupling case: Plugging this ansatz inside Eq. (5.39) and differentiating twice with respect to χ, we find with ζ = χ/Y α . To ease the counting of Y powers, we have expanded the right hand side, equal toᾱ s (ζY 1/2+α + ρ s )Y α G for large Y .
The strategy is to find the leading power in Y in Eq. (5.51) in order to get a simpler second order differential equation for G. The r.h.s. being proportional toᾱ s , it must contribute to this differential equation, and therefore it must contain this leading power. Since the first term b 0 Y α−1 G cancels against a similar term in the l.h.s., the leading power is associated with the second term, proportional to Y 2α−3/2 . As we aim at finding a second order differential equation, the leading power in the coefficient of G must also be proportional to Y 2α−3/2 . This homogeneity condition 2α − 3/2 = −1 − α yields α = 1/6, implying that the leading power in Eq. (5.51) is Y −7/6 . As noted in the previous section, the diffusion power α = 1/6 is not the same as the 1/2 found in the fixed coupling evolution, showing that in the running coupling case, the evolution equation does not belong to the same universality class as the FKPP equation.
For consistency, one has to include a term of order Y −5/6 in the development ofρ s , namelyσ with an unknown constant δ 1 to be determined. Indeed, such term generates contributions of the same order as Y −7/6 . Gathering all terms proportional to Y −7/6 provides the differential equation Notice that the ζ dependence of the right hand side of this equation comes from the ρ argument in the QCD coupling. If we had replacedᾱ s (ρ) byᾱ s (ρ s (Y )) in the evolution equation, this term would not be there, changing the value of the constant δ 1 that we now determine.
The only solution of the Airy-type differential equation (5.53) which satisfies G(ζ) = √ b 0 ζ + O(ζ 2 ) in order to recover the leading behaviour of the scaling solution f near the interior of the front and with vanishing boundary conditions at ζ = +∞ is The constant δ 1 is then fixed by these boundaries conditions: The study of the diffusion around the scaling limit in the leading edge domain enables us to recover the asymptotic expansion of ρ s (Y ): and provides the first correction G to the scaling limit. Compared to the fixed coupling case, one notices that the first two terms in the development of ρ s does not depend on the linearization of the evolution equation.

Corrections in the interior and on the edge of the wavefront
So far, we have essentially recovered the results obtained in section 5.1 for a linear saturation boundary. The main advantage of the mathematical approach of wave front formation is that it provides a systematic way of calculating the higher orders in the asymptotic development of ρ s (Y ) and the sub-asymptotic corrections ofq(Y, ρ) even in the presence of a non-linear saturation boundary.
Front interior expansion. The series expansion of the diffusive solution (5.54) is instructive. One finds This expansion enables one to infer the form of the front interior expansion: at fixed χ, the term in ζ 3 brings a power Y −1/3 , the term in ζ 4 brings a power Y −1/2 and so on. We therefore define the front interior expansion aŝ The rescaling Y → 4b 0 Y and χ → √ b 0 χ is purely conventional and simplifies the expressions of the functions f n . The function f 0 has already been computed. From the expression of the scaling limit f (χ), one gets f 0 (X) = 1 + X.
Using the evolution equation and the following development ofρ ṡ where the presence of the corrections in 1/Y , 1/Y 7/6 and 1/Y 4/3 will be justified a posteriori from our calculation of the next terms in the leading edge expansion, we find that the functions f n follow an infinite hierarchy of second order differential equation of the form f n (X) = .... As in the fixed coupling calculation, this hierarchy can be solved iteratively since the system is triangular (the right hand side depends only on f i (X) with i < n). Concerning the initial conditions, the definition of the saturation boundary yields f n (0) = 0 for all n ≥ 1. The conditions on the first derivative is obtained thanks to the differential equation in integral form. The first two terms read One observes that the leading power in the polynomial functions f 0 , f 1 and f 2 are included in the series expansion of G displayed in Eq. (5.57). In the front interior expansion (5.58), the first six terms (up to n = 5) are universal. By universal, we mean that they do not depend on the initial conditions for the evolution nor on the constant term in the development of ρ s (Y ) at large Y . All the functions f n (X) for 0 ≤ n ≤ 5 are provided in appendix A.
Leading edge expansion. The front interior expansion does not enable to fix the value of the coefficients δ i in the asymptotic expansion ofρ s . These coefficients are determined by matching the leading edge expansion with the front interior one. In the running coupling evolution, the leading edge expansion takes the form . After a rather tedious calculation similar to the one leading to the solution G −1 , one gets the following differential equation for G 0 : by looking at the coefficient in front of the Y −4/3 power on both side of the equation. The right hand side depends on the coefficient δ 2 in front of the Y −1 power in the asymptotic development of ρ s , it is then a consistency requirement to include such a term, and justify a posteriori the form of (5.59).
Contrary to the calculation of section 4.2, the homogeneous equation is the same as the one satisfied by G. Since the inhomogeneous right hand side is known, this equation can be solved by the method of variational parameters. The initial conditions for G 0 are G 0 (0) = 1 and G 0 (0) = 0 in order to match with the front interior expansion which has no term of order χY −1/6 . These two conditions fix the two constants of integration. The constant δ 2 is then determined by demanding the solution G 0 to decay exponentially at large ζ. We find that the coefficient δ 2 in front of the ln(Y ) term reads The coefficient in front of the ln(Y ) term is different from the coefficient 1/4 found in [23] by numerically solving the linearized evolution equation forq 2 . As in the fixed coupling case, the non-linearity of the saturation boundary brings a sizeable correction in front the logarithmic term in the asymptotic development of the saturation scale. Parametrically, this correction is of order b 0 ln(Y ) ∼ α s (ρ s )Y ln(Y ) and is then larger than pure single log corrections, or order α s Y . The sub-asymptotic correction G 0 (ζ) can be expressed in terms of the Airy function and its derivative, Ai (s(ζ)) + − 7 12 Ai(s(ζ)) , It is quite remarkable that the dependence upon the QCD constant b 0 enters only through this shift function s(ζ). This feature does not persist in the higher orders of the leading edge development. One can also verify that this function admits the series expansion and therefore accounts for the sub-leading powers of the front interior expansion. 2 We have also checked that our analytic approach enables to recover this 1/4 coefficient for the linearized evolution equation, cf. section 5.3.
Universality. Finally, one may wonder which terms in the asymptotic expansion of are universal. By universal, we mean independent of the initial condition forq. In particular, the integration constant κ is not determined by the leading edge regime and depends on the initial condition. Therefore, the terms in the development of ρ s orq(Y, ρ) at large Y which depends on κ are not universal. Expanding α s (ρ = ζY 2/3 + ρ s ) in powers of Y , one gets When this development is multiplied by the leading edge expansion, the smallest power of Y which involves the coefficient κ is Y −2+1/6 = Y −11/6 . This is the power which determines the differential equation satisfied by G 3 . Therefore, the function G 3 is not universal while G 1 and G 2 are. These functions can be determined by following the same method as in the computation of G 0 . The coefficients δ 3 and δ 4 in the development (5.59) ofρ s are obtained from the matching of G 1 and G 2 with the front interior expansion and from the boundary conditions in ζ = ∞: The two functions G 1 and G 2 are provided in appendix A. Last, one notices the presence of a term of order ln(Y )/Y 2 in (5.68). In principle, such term would spoil the shape of the leading edge development, introducing a contribution of the form Y −1/2 ln(Y )G 3 (ζ). Yet, it is not the case 3 as one can prove that the only solutionG 3 consistent with the front interior expansion isG 3 = 0, provided that the asymptotic expansion ofρ s readṡ and is indeed not universal since it is κ and ρ 0 dependent. For ρ s , our final result including all universal terms (except for the unknown integration constant κ) in the large Y expansion is One notices that the b 0 dependence cannot be absorbed into a redefinition of the variable Y → b 0 Y because of the b 0 dependent terms inside the parenthesis of the ln(Y ) and Y 1/3 terms.

Comparison with the linearized evolution equation
We now discuss the asymptotic expansion ofq(Y, ρ) and ρ s (Y ) for the linearized evolution equation. In this scenario, one can also apply the analytic techniques of front propagation into unstable states to derive the large Y development of these quantities. The only differences are the change of scaling variable 73) and the definition of the saturation scale which is not defined by an implicit relation but rather as One can then compute the front interior, the leading edge and the ρ s asymptotic expansions as in the non-linear case. The formula for the front interior and leading edge functions are provided in appendix A. They display interesting scaling properties in terms of the variable b 0 Y and s respectively. This feature persists for the asymptotic expansion of ρ s which reads Thus, the presence of the terms linear in b 0 in the coefficient of the asymptotic series in Eq. (5.72) is a consequence of the back-reaction of the quantum evolution ofq on the saturation boundary. An important remark concerns the absence of term of order ln(Y )/ √ Y in this development. The reason is that for the linear evolution, the right hand side of the differential equation satisfied by the functions of the leading edge development is given by the expan- In other words, the scale ρ s (Y ) is simply replaced by Y due to the absence of back-reaction, so that there is no spurious ln(Y ) terms in the power development above. For the same reason, the non-universal contribution of order Y 1/2 in Eq. (5.75) does not depend on κ and simply reads In spite of being non-universal, this contribution enables to accurately describe the function ρ s (Y ) down to small Y values (Y ∼ 1 ÷ 2).

Numerical study: convergence of the asymptotic espansion
The expansion (5.72) and its truncations up to the orders O(Y −1/2 ) and O(Y −1 ) are compared to numerical computation of ρ s (Y ) in Fig. 8b. Contrary to the fixed coupling case, the convergence of the this asymptotic expansion at moderate values of Y is much slower. This is essentially a consequence of the weaker power Y −1/6 of the scaling deviation compared to the fixed coupling case Y −1/2 , due to the fact that the two evolution equations do not belong to the same universality class. Because of the large b 0 dependent contribution in front of the ln(Y ) and Y −1/3 terms, the coefficients of the power expansion in Y ofρ s behave like the coefficients of an asymptotic series, and therefore there is no convergence at small or moderate values of Y .
The convergence of the leading edge expansion is shown in Fig. 8a. The red curve includes the first two terms G −1 and G 0 , while the green curve includes all the universal terms in both the leading and front interior expansion. For Y = 100, a rapid convergence of the series is observed. Unfortunately, for the same reason as for ρ s (Y ), at low and moderate values of Y , the truncated series has a pathological behavior due to both the slow convergence of the series and the linear terms in b 0 which behave like a divergent asymptotic series.
In the next section we will address this short coming at moderate values of Y . We abandon, in particular, the leading edge expansion for a Taylor series in the vicinity of the saturation line, x = 0 , i.e., ρ ∼ ρ s (Y ).

Phenomenology of transverse momentum broadening
In this section, we study the transverse momentum distribution after resummation of the leading radiative corrections, including running coupling effects. We aim at providing analytic expressions for the k ⊥ distribution, that may be used for the phenomenology of transverse momentum broadening in heavy-ion collisions or as initial conditions for the small-x evolution of gluon distributions in large nuclei.

Taylor expansion about ρ s
Although the leading edge expansion forq is essential to determine systematically the sub-asymptotic corrections to ρ s , we have shown that for moderate values of Y it fails to converge, and therefore cannot be used at small values of Y = 1 ÷ 5, the typical values relevant for phenomenology. This limitation is only apparent. Indeed, the leading edge expansion applies to large x values away from the saturation regime. Hence, if one is mostly interested in finite region near the saturation scale then a Taylor expansion around x = 0 should be enough to achieve a desirable accuracy: with x = ρ − ρ s (Y ) as usual. This is not to say that the leading edge expansion is not useful, quite the contrary. Indeed, we will show that the coefficients of this Taylor expansion depends only on the function ρ s (Y ) which itself is determined by the leading edge expansion. Therefore, if we have a good analytic knowledge of this function at small Y , we may expect that the resulting TMB distribution will be close to the exact numerical result. It turns out that at small Y , the effect of the non-linearities in the evolution equation are mild, so that the function ρ s (Y ) in the non-linear case is in fact close to ρ s (Y ) in the linear case. Since we have a very good analytic control of ρ s (Y ) in the linear case, we will use its form as our analytic input inside the Taylor expansion that we now detail. As reported in Sec. 2.2, when using the functionsq(Y = ln(L/τ 0 ), ρ) and ρ s (Y ) to evaluate the transverse momentum distribution, one must distinguish between the dense regime ρ < ρ s (x < 0) and the dilute one ρ > ρ s (x > 0). First, we want to determine the Taylor expansion of F(x, Y ) = F > (x, Y ), to the right of the saturation line, up to second order in x around x = 0 (ρ = ρ s (Y )). As a matter of fact, there is a one-to-one correspondence between ρ s (Y ) and the Taylor coefficients of F > (x, Y ) at x = 0. Indeed, plugging the definition of F > inside the differential equation satisfied byq(Y, ρ) one gets Using the fact that F > (0, Y ) = 1, from which we also deduce thatḞ(0, Y ) = 0, Eq. (6.3) is readily solved yielding To obtain the second derivative we need to differentiate Eq. (6.3) w.r.t. x Evaluating the latter at x = 0 and using Eq. (6.5), we obtain As a result we obtain with ρ s and its derivatives evaluated at Y = ln(L/τ 0 ). Higher order derivatives, F Let us turn now toq(ρ) =q < (ρ) to the left of the saturation line. By definition we haveq , and Taylor expanding with respect to x, one finds that In the end, one can approximate the functionq < usinĝ If one expands again this result in powers of x, one notices that the coefficients of x 0 and x are equal to those ofq > , but not the coefficient of x 2 . It means that the functionq(ρ) we will employ in our analytic results is continuous and derivable, but not twice differentiable.

Numerical results
The formula forq(Y, ρ) obtained from the Taylor expansion approach depends on ρ s (Y ) only. If one aims at quantifying the effects of the quantum corrections for realistic values of Y , say Y = 2 ÷ 5, one needs an analytic expression for ρ s (Y ) that correctly describes this range of values. In Fig. 9a, we observe that at moderate Y , the non-linear (blue curve) and the linear (red curve) evolution of the saturation scale are very close. On the other hand, even though the effects of the non-linearities are mild, we are already in the universal regime which is not driven by the initial condition. One can then take advantage of this fact, since contrary to the non-linear evolution, the asymptotic expansion of ρ s (Y ) is convergent in the linear case, even at small values of Y . The convergence of the asymptotic development (5.75) is demonstrated in Fig. 9b, where the dashed curve labeled O(Y −1/2 ) also includes the term given by Eq. (5.77). In these analytic curves, the value of the unknown integration constant κ is determined by fitting the large Y tail of the numerical data. In the non-linear case, as we have shown in the previous section, the development is divergent and does not describe the numerical data for Y smaller than 10. This is clear from the grey curve in Fig.9a, when compared to the numerical result in blue.
In the linear case, the asymptotic expansion of ρ s (Y ) at O(Y −1/2 ) accuracy provides a very good approximation of both the linearized and "exact" ρ s (Y ): in Fig. 9, the dotted black curve overlaps with the red curve even at Y of order 1 and is also close to the blue curve. The mathematical reason is that the linear development at large Y is related to the non-linear one by dropping the problematic terms (proportional to b 0 , induced by the back-reaction of the quantum evolution to the saturation boundary) which makes the series divergent as Y becomes smaller.
Since at moderate Y values, the effect of the non-linearities on the saturation scale ρ s (Y ) are mild, it is legitimate to use the expressions (5.75) and (5.77) as our analytic input for ρ s in the Taylor expansions given by Eq. (6.7)-(6.10). We have now all the ingredients to compare the functionq(Y, ρ) obtained by solving the non-linear evolution equation numerically with our analytic expressions (6.7)-(6.10). This is shown Fig. 10a: the red curve is our numerical result forq(Y, ρ) as a function of x = ρ − ρ s , while the dashed black curve correspond to Eq. (6.7)-(6.10). This analytic approach is very conclusive and can be systematically improved by including more terms in the Taylor expansion in order to describe the large x domain.
In Fig. 10a, the grey curve is the numerical result forq(Y, ρ) from the linearized evolution equation. In contrast with what we observed for ρ s (Y ), we notice that the non-linearity has an important effect onq: it slows down the evolution. The difference comes from the coefficients of the Taylor expansion that we have established in the previous sub-section: for the linear evolution equation, the dependence upon ρ s of these coefficients is not the same. In particular, the first derivative in x = 0 isρ s − 1 which is significantly larger than 1 − 1/ρ s at moderate values of Y . Hence, even though we observed that the asymptotic development of the linearized ρ s is a good proxy for the "exact" ρ s for Y of order 1, the linearization turns out to be a bad approximation for the functionq(Y, ρ) and consequently, for the TMB distribution itself.
The resulting TMB distribution is shown in Fig. 10b for two values of Y , Y = 2 and 4. The dashed curves correspond to the analytic expressions after Fourier transform of the dipole S-matrix. For Y = 4, the agreement is excellent, and even for Y = 2, our formulas correctly captures the general trend of the distribution. We point out that the oscillatory behaviour at large k ⊥ is a consequence of the discontinuity of the second derivative of q(L, k 2 ⊥ ) with respect to k 2 ⊥ in Q 2 s .

Physics discussion
We now detail the physical interpretation of the asymptotic limit of the quenching parameter when the system size L goes to infinity. We rely mainly on the equation established in section 4.2.1. The discussion will be divided into two subsections: one related to the behaviour of the asymptotic TMB distribution around its peak, and an other about its large k ⊥ tail.

Anomalous diffusion and Lévy flights
Using Eq. (4.23) and ρ s (Y ) = cY at large Y , one can obtain the scaling limit ofq in the dense regime that controls the peak of the k ⊥ distribution. Indeed, we havê For k 2 ⊥ ≤ Q 2 s (L), we need to evaluate Y along the saturation line as shown by Eq.(2.26). Slightly abusing the notation, since we now name Y = ln(L/τ 0 ), we havê where we have used f (0) = 1 and Y s (ρ) ≈ ρ/c in the scaling limit (we recall that β = (c − 1)/(2c) and c = 1 + 2 ᾱ s +ᾱ 2 s + 2ᾱ s is the velocity of the traveling wave front). Plugging this expression inside the dipole S-matrix and transverse momentum distribution, one obtains the two dimensional Fourier transform of a "stretched exponential": Therefore, the leading effect of the radiative corrections for large system size is a modification of the power exponent of the effective dipole size |x ⊥ |. In the running coupling case, one would obtain similar results, with the value of β ≈ √ᾱ On top of this anomalous scaling of the "moments" of the k ⊥ distribution, we emphasize that its whole shape around the peak differs significantly from the tree-level result. This also applies even at moderate values of the system size and not only in the rather formal asymptotic L → ∞ limit, as shown in Fig. 4, where one observes that the resummed distribution is significantly broader than the tree-level one.

Large-k ⊥ tail of the TMB distribution
In the previous section we have discussed the behavior of the distribution around its peak that is controlled by the saturation scale Q s (L). At large k ⊥ , away from the saturation line where the physics is dominated by a single hard scattering there is an interesting relation between the large k ⊥ tail of the TMB distribution and the medium gluon distribution associated with the correspondence (8.1). To extract the large k ⊥ behavior of the TMB distribution, we expand the exponential in the definition of the dipole S-matrix In the limit k T → +∞, the z integral is dominated by z ∼ 1, and we thus have ln 1/z 2 ln k 2 ⊥ . It allows us to expandq aŝ we obtain the formula [19] 11) where xG N is the gluon distribution function defined in Sec. 8.2 and x = max(τ 0 /L, k 2 ⊥ /2P + T ) (the value of x used in this expression will be commented in the next section). In the second line, we have assumed that k 2 ⊥ dependence of dq/d ln(k 2 ⊥ ) is relatively slow compared to the power law behavior due to the 1/k 2 ⊥ factor. To measure the qualitative effect of the k 2 ⊥ dependence ofq as a result of quantum corrections, it is convenient to parametrize the deviation to the Rutherford behavior as follows: When E ω c , one must distinguish several different regimes at large k ⊥ , as shown Figure 11. In the domain Q s k ⊥ Q 2 s /µ, we are in the extended geometric scaling window, whose behaviour is driven by the asymptotic limit (Y → ∞) of the quenching parameterq(Y, ρ). This asymptotic limit satisfies geometric scaling, in so far asq(Y, ρ) is a function of x = ρ − ρ s = ln(k 2 ⊥ /Q 2 s ) only. Indeed using Eq. (8.12), one finds that for In the extended geometric scaling window, the TMB distribution behaves like a power law, with slower decay than the Rutherford one. This heavy tail is also characteristic of Levy flight processes. At fixed coupling at least, the analogy with Lévy flights enables to understand both the anomalous behavior of the diffusion process, since Q s ∝ L 1 2 + √ᾱ s instead of L 1/2 and the heavy-tailed distribution. One could argue that even at tree-level, the TMB distribution displays a heavy tail with a Rutherford like decay. However, we emphasize that the heavy tail highlighted in this paper is entirely a consequence of the radiative corrections. Indeed, even in the case of a constant initial conditionq (0) (ρ) =q 0 (the so-called "harmonic approximation") for which the tree-level distribution is gaussian and therefore decays exponentially at large k ⊥ , the resummation of self-similar gluon fluctuations leads to a heavy, power law tail given by Eq. (7.14).
Beyond the extended geometric scaling window, for k ⊥ Q 2 s /µ, but for k ⊥ ≤ E/L, the k ⊥ -distribution becomes sensitive to the double log of the DGLAP evolution, and the exponent D(k 2 ⊥ ) reads The function D is therefore slowly tending to 0 at large k ⊥ , meaning that one recovers asymptotically the Rutherford 1/k 4 ⊥ power law, as illustrated Fig. 11 between the scales Q 2 s /µ and E/L.

Energy dependence of the quenching parameter
So far, we have considered the propagation of a highly energetic parton, with E = P + ω c . In this section, we comment on the opposite regime E ω c . In this limit, the characteristic quantum diffusion time 2E/k 2 ⊥ of the incoming parton wave-function can become smaller than the system size L even in the saturation regime where k 2 ⊥ ∼ Q 2 s ≈q 0 L. In order to understand what the dominant radiative corrections are in this case, we make a short detour by the operator definition ofq and its link with parton distribution functions (PDFs). Then, we address the effect of the leading radiative corrections when E ω c , and derive the typical energy dependence of the quenching parameter in this regime.

8.1
On the x-dependence ofq q and gluon PDF. The quenching parameter is to some extend related to the gluon parton distribution function. However, the connection is not straightforward and certainly the two quantities are not equivalent across the full k ⊥ spectrum. To begin with, PDF's are defined in the dilute regime of weakly interacting partons while the quenching parameter q is sensitive to non-linear or saturation effects. Bearing in mind these differences, it is instructive to explore the identification for cold nuclear matter  Figure 11: An illustration of the various regimes of the transverse momentum broadening distribution P (k ⊥ ). In the non-linear regime k ⊥ Qs = √q L the distribution is dominated by multiple soft scattering and the distribution exhibits geometric scaling, i.e., it is only a function of the scaling variable k ⊥ /Q s . The regime Q s k ⊥ Q 2 s /µ corresponds the extended geometric scaling window that is characterized by a heavy tailed distribution. For k ⊥ > Q 2 s /µ the dynamics is linear and thus, not sensitive to saturation physics. k ⊥ = E/L scale marks the transition from k − ≡ xP − ∼ 1/L to k − ≡ xP − ∼ k 2 ⊥ /2E. The latter relates to the standard Feynman x of the medium PDF.
The function xG N (x, Q 2 ) is the the gluon PDF of a nucleon target, which reads in the parton model [69] xG is the field strength tensor and [x + , 0 + ] is a gauge link connecting the points (x + , 0 − , 0) and (0 + , 0 − , 0), in light-cone variables. In doing so we recover the Glauber-Mueller formula for the forward dipole scattering amplitude [16,70]. However, there is an ambiguity in the Feynman x-dependence of the gluon distribution in the idenfication (8.1). The latter can be resolved by computing quantum corrections. At sufficiently low k ⊥ the one loop correction is dominated by soft and collinear double logarithmᾱ where τ ∼ 1/k − is the formation time of the soft gluon radiation and k ⊥ its transverse momentum. Including these radiative corrections through the renormalization of the gluon PDF makes (8.2) Q 2 -dependent. Similarly,q acquires a Q 2 dependence through the resummation of double logarithmic radiative corrections like (8.3). In the case ofq, this Q 2 dependence should not be confused with the "tree level" Q 2 -dependence coming from the upper limit of the k 2 ⊥ integration whenq is defined as k 2 ⊥ /L, as in [47,70]. In that case, the Q 2 dependence is a higher twist effect resulting from the Coulomb tail of the single gluon exchange in the t channel which makes the second moment of the TMB distribution divergent.
In a similar fashion, the x-dependence ofq can be understood from the τ max dependence induced by higher order corrections. To set the typical value of τ max , one must distinguish two regimes. Indeed, from the on-shellness condition of the radiated gluon, one gets the relation k − = k 2 ⊥ /(2k + ). Since k + < P + , the latter equality gives the kinematic constraint τ max ∼ 2P + /k 2 ⊥ . 4 In the TMB process, the two regimes are then either τ max > L or τ max < L.
In the case τ max > L, the largest formation time is actually set by the medium size L, namely 1/k − ∼ 1/(xP − ) < L, and one recovers after resummation the double logarithmic limit (2.31) (with Y = ln(L/τ 0 )) of the evolution considered in this paper. Indeed, fluctuations larger than L are strongly suppressed due to the LPM effect owing to the fact that these long-lived gluons do not resolve the medium from the hard scattering that originates the jet. Notice that this constraint does not apply in the case of an asymptotic quark scattering off a shock-wave where the so called small-x gluons can stretch beyond the extent of the target up to τ max ∼ 2P + /k 2 ⊥ > L. Effectively, the relation τ max ∼ L translates into the upper bound L for the x + integration in the operator definition ofq, On the other hand, when τ max ∼ 2P + /k 2 ⊥ L the gluon fluctuation is not sensitive to the size of the system and we must recover the standard gluon PDF where xP − ≡ k 2 ⊥ /2P + and P − ∼ T . The small and large x regimes ofq that are encompassed by Eq. (8.4) and Eq. (8.5), respectively, can be combined in This definition differs from other definitions encountered in the literature [37,72,73], in whichq is defined as the second moment of the TMB distribution at leading twist. We believe that our calculation is more natural given thatq appears in the unintegrated distribution and thus depends locally on k ⊥ . To sum up this discussion, our main findings are that the x-dependence, or equivalently the τ -dependence of the quenching parameter should be given by or in terms of x Differences between the evolution ofq and DGLAP/BFKL. Based on the above insight on the relation betweenq and the gluon PDF, let us now analyze more closely the double logarithmic structure in Eq. (8.3) in order to identify the major differences with the two main regimes of QCD evolution, namely, DGLAP and BFKL. First, DGLAP evolution equations resum powers of the single logarithm ln Q 2 /µ 2 which in our context would be ln k 2 ⊥ /µ 2 , while assuming that x is of order 1 such that there is no need for resumming ln 1/x powers. This criterion is obviously met at large enough k ⊥ at the end of the power tail as shown in section 7.2 (cf. also figure 11 and the shape of the k ⊥ distribution between Q 2 s /µ and E/L). Now when x 1 in Eq. (8.6) but at the same time we have xP − L 1, i.e, k ⊥ > P + /L, we would have at next to leading order After resummation, usingq one finds that the deviation D to the Rutherford behaviour defined in Sec. 7.2 is given by D(k 2 ⊥ ) 4 (ρ E − ρ)/ρ. This is illustrated in figure 11 in the domain k ⊥ > E/L (when k ⊥ gets close to √ ET , one approaches the kinematical limit beyond which the eikonal approximation is no longer valid).
For k ⊥ < P + /L the quantum phase in Eq. (8.6) is no longer relevant and the integral over x + must be cut-off at L. This yieldŝ In this regime the collinear logarithm ln k 2 ⊥ µ 2 is always larger that the soft logarithm ln L τ 0 and we may wonder whether the relative importance of these logs would be reversed, in which case the Regge kinematics where single logs of 1/x are resummed would be more appropriate. However, saturation effects are relevant precisely in the regime where both logs are equally important. Indeed, multiple scattering responsible for the unitarization of the dipole amplitude are effective when k 2 ⊥ = Q 2 s =qL. As a result, we obtain the double log structureᾱ s ln 2 L τ 0 . Unless we are interested in the deep saturation regime k ⊥ Q s a single log approachà la BFKL is not necessary.
On the other hand, in the case of DIS in the standard Regge kinematics the saturation scale scales as Q 2 s (x) ∼ x λ where the anomalous dimension λ ∝ α s 1 [74], we thus have ln Q 2 s (x) = α s ln 1/x ln 1/x, where we readily see that the collinear log is suppressed compared to the soft log by the coupling constant.
In contrast, in the case of momentum broadening the leading order comes with a factor Q 2 s ∼ L ∼ 1/x (where we used that xP − ∼ 1/L) and quantum evolution will only slightly depart from this behavior. As we have seen, we have in particular Q 2 s ∼ L 1+2 √ᾱ s at leading logarithmic accuracy. In other words, the problem of transverse momentum broadening near the saturation line is of double logarithmic nature and hence, does not favor neither DGLAP nor BFKL evolution equations. This justifies a posteriori why the dilute regime is sufficient for the qualitative analysis of the phase space.

E
qL 2 and energy dependence of the saturation scale We now return to discussing the qualitative consequences of the relation (8.7) on the energy dependence of Q s . In DLA, the asymptotic behavior ofq readŝ as can be shown from Eq. (2.29) 5 . In this expression, the value of τ is bounded by the minimum between the jet path length L and the quantum diffusion time given by 2E/k 2 ⊥ . In the high energy limit, i.e. P + = E ω c ≡qL 2 , which is the main focus in this paper there is a potentially large phase space in which min L, 2E This phase space is given by Q 2 s ≤ k 2 ⊥ ≤ 2E/L. In this regime, using the defining equation for Q s and Eq. (8.12), one readily finds ρ s (Y ) = c Y , c = 1 + 2 ᾱ s +ᾱ 2 s + 2ᾱ s , (8.14) corresponding to the superdiffusive scaling.
Let us now explore the opposite limit, which corresponds to E ω c . In this case we have 1/τ = k 2 ⊥ /2E 1/L not only at high k ⊥ but also deep inside the saturation regime. We can determine the behavior of the saturation momentum from its defining equation and the asymptotic form (8.12): Defining the variable ρ E = ln(2E/(q 0 τ 2 0 )) ∼ ln 2ET /µ 2 , the above equation can be rewritten as Solving the quadratic equation for ρ s and retaining the larger solution we obtain (8.17) 5 We have dropped sub-leading power prefactors. In the limitᾱ s 1 andᾱ s ρ E Y , we may simplify further that is This is the typical leading order result, which does not display the anomalous behavious. As a function of L for a fixed parton energy E, the saturation momentum recovers its linear scaling with the system size for L large enough, L 2E/q 0 . The leading behaviour of the saturation scale as a function of the system size Y = ln(L/τ 0 ) is shown in Fig. 12a, where one observed the two regimes (8.14) and (8.17) with a transition around the scale 2E/q 0 . Conversely, for a fixed system size L, Q s exhibits an energy dependence below the scale ω c . This dependence, given by Eqs. (8.14)-(8.17), is represented on Fig. 12b where ρ s is plotted as a function of ln(2E/(q 0 τ 0 )) ∼ ln(E/T ).

Summary
In this article, we investigate quantum corrections to transverse momentum broadening in the double-logarithmic approximation by solving analytically and numerically the nonlinear evolution equation for the quenching parameterq recently put forward.
The effects of quantum evolution are three-fold: i) They cause a heavy tailed distribution, akin to Lèvy random walks, to form at large Y = ln L/τ 0 . i) The TMB distribution loses sensitivity to initial condition given by the tree-level and tends to a universal distribution that can be computed analytically. iii) The TMB distribution obeys a geometric scaling in an extended region above the saturation scale, i.e., Q s k ⊥ Q 2 s /m D .
Our analytic approach is based on an asymptotic expansion of the solution for large Y in both fixed and running coupling using techniques that we borrowed from traveling waves analyses and gluon saturation physics. In the regime where the fast parton energy tends to infinity, we show that the transverse momentum distribution quickly reaches a universal regimes where it is approximately a function of a single scaling variable k ⊥ /Q s (L) with mild scaling violations that can be computed systematically. We have in particular computed the first six terms in the asymptotic expansion of Q s (Y ).
Furthermore, we have derived new results for the running coupling case that differs substantially from the fixed coupling case discussed in a previous work. We have in particular provided a formal proof for the asymptotic formula for Q s (Y ) conjectured in Ref. [23] and computed the corrections to their results caused by the non-linearity of the saturation boundary. The running coupling exhibits a weak scaling which can be understood approximately from the fixed coupling analysis as an additional Y dependence through the coupling constant that enters the exponent of the power spectrum.
Our approach consists in searching for diffusion-like solutions near the exact scaling solutions away from saturation line. This region at the edge of the wave front leads the propagation of the front and in turn, is responsible for setting the overall speed of the wave given byρ(Y ).
In the fixed coupling case, where geometric scaling is asymptotically exact, we were able to obtain a very good agreement between our universal asymptotic expansion with the numerical computation for Q s (Y ) down to small values of Y ∼ 2 − 5. On the other hand, the asymptotic expansion in the running coupling case, when the non-linear effects are accounted for, does not converge at low values of Y . We nevertheless observed that the expansion in the linearized equation provides a good approximation.
The phenomenological interest of the analytic formulas derived in this paper is twofold. First, they can be used to estimate at low numerical cost the TMB of an energetic jet propagating through the quark-gluon plasma formed in heavy-ion collisions. In the case of the dijet azimuthal asymmetry, computed in the Sudakov formalism in [13,14], TMB enters in the calculation as a single parameter, the average transverse momentum squared k 2 ⊥ , or the saturation scale Q s , which is fitted to experimental data. The relation between this parameter and the medium physical properties can then be obtained from our expressions for Q s (L). It would also be interesting to experimentally observe the super-diffusive regime by scanning the saturation scale over various system sizes (e.g. various centrality classes or nuclei).
Another possible application of our results pertains to QCD at small-x. At high energy and fixed Q 2 , the building block of the fully inclusive DIS cross-section is the dipole S-matrix (2.2) whose evolution with Bjorken-x is governed by the non-linear BK equation. In phenomenological studies, the BK equation is usually solved using the McLerran-Venugopalan model [75,76] as the initial condition at moderate values of x = x 0 ∼ 0.01. The MV model is analogous to the tree-level form of S(x ⊥ , L) that gives the tree-level TMB distribution. This paper shows that including the leading radiative corrections enhanced by the nucleus size, of order α n s ln 2n (A 1/3 ) for all n ≥ 0 leads to a universal distribution which becomes insensitive to the detail of the tree-level physics. This motivates the use of our analytic expression as a new initial condition for the BK equation. This approach would differ from other attempts to go beyond the MV model such as [77,78] which focuses on power of 1/A suppressed corrections in the MV effective action or [79,80] which considers only a single gluon emission from a valence quark in the computation of the two-point correlator.
We conclude this summary with a brief perspective on future studies. The asymptotic expansions calculated in this paper are valid in the double logarithmic limit. We expect the single logarithmic corrections to change some of the coefficients in this series, as shown in [60] in the context of the BK or BFKL evolution. For the jet quenching parameter problem, the resummation of the single logarithmic corrections raises additional questions related to the proper kernel to be used in the evolution or the exponentiation of the single log terms. Nevertheless, it is crucial for precision phenomenology to go beyond the present results. In order to draw a complete picture of the transverse momentum broadening in a dense QCD medium over all transverse momenta, such single logarithmic resummation should be matched with NLO results with exact kinematics [20], NLO corrections to the collision kernel [81][82][83] and a proper determination of the non-perturbative small k ⊥ domain (e.g from lattice calculations [84,85]). Regarding the phenomenological applications to heavy-ion collisions, one could also investigate the relative importance between higher order corrections and the effects of inhomogeneities in the plasma as computed in [86].
with χ = (ρ − ρ s (Y ))/Y 1/2 . Only the first five terms in this sum are universal: f 0 (X) = 1 + X , (A.2) One notices that all these f n functions are functions of X = √ b 0 χ and b 0 , f n → f n (X, b 0 ). One can show that in the linear case, the functions f n are given by setting b 0 = 0 in the second argument, i.e. f n = f n (X, 0) (the variable χ becomes also χ = (ρ − Y )/ √ Y ). Therefore, the b 0 scaling associated with Eq. (A.1) is exact for the linear evolution equation.

The universal leading edge functions
The leading edge expansion resums all the leading powers in the front interior expansion, as explained in the main text. In mathematical terms, the leading edge series readŝ The functions P n and Q n are functions of the variable s and b 0 through the coefficients of these polynomial functions. As for the front interior expansion, the functions P n and Q n for the linearized evolution equation are obtained by setting b 0 = 0 in all these coefficients.