# Transverse momentum broadening from NLL BFKL to all orders in pQCD

Paul Caucal, Yacine Mehtar-Tani

► To cite this version:

Paul Caucal, Yacine Mehtar-Tani. Transverse momentum broadening from NLL BFKL to all orders in pQCD. Physical Review D, 2023, 108 (1), pp.014008. 10.1103/PhysRevD.108.014008. hal-03793084

HAL Id: hal-03793084

https://hal.science/hal-03793084v1

Submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Transverse momentum broadening in large media
from NLL BFKL to all orders in pQCD

Paul Caucal<sup>1,2,*</sup> and Yacine Mehtar-Tani<sup>1,3,†</sup>

¹Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA

²SUBATECH UMR 6457 (IMT Atlantique, Université de Nantes, IN2P3/CNRS),
4 rue Alfred Kastler, 44307 Nantes, France

³RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973, USA

<figure><img src="image_2.png" /></figure>

(Received 22 October 2022; accepted 8 June 2023; published 11 July 2023)

We study, to all orders in perturbative QCD, the universal behavior of the saturation momentum $Q_s(L)$ controlling the transverse momentum distribution of a fast parton propagating through a dense QCD medium with large size $L$. Due to the double logarithmic nature of the quantum evolution of the saturation momentum, its large $L$ asymptotics is obtained by slightly departing from the double logarithmic limit of either next-to-leading log (NLL) BFKL or leading order DGLAP evolution equations. At fixed coupling, or in conformal $\mathcal{N} = 4$ SYM theory, we derive the large $L$ expansion of $Q_s(L)$ up to order $\alpha_s^{3/2}$. In QCD with massless quarks, where conformal symmetry is broken by the running of the strong coupling constant, the 1-loop QCD $\beta$-function fully accounts for the universal terms in the $Q_s(L)$ expansion. Therefore, the universal coefficients of this series are known exactly to all orders in $\alpha_s$.

DOI: 10.1103/PhysRevD.108.014008

# I. INTRODUCTION

The study of the suppression and modification of jets in heavy-ion collisions at RHIC and LHC [1–9], commonly referred to as “jet quenching” [10–12], aims at probing the quark-gluon plasma at various scales as well as nonequilibrium dynamics of QCD. This physics will play a central role in the sPHENIX program at RHIC [13] and the upcoming run 4 at the LHC [14]. In that context, the phenomenon of transverse momentum broadening (TMB) of jets in the QGP is of major importance as it controls many jet quenching related observables measured in heavy-ion collisions. For instance, TMB is responsible for the jet energy loss by deflecting soft medium-induced gluons at larger angles than the jet cone size $R$ [15] leading to the suppression of the jet cross section in nucleus-nucleus collisions [16–19]. Another important and historical signature of jet quenching is the dijet azimuthal asymmetry [20–23]. This observable is believed to be sensitive to TMB of jets propagating in the quark-gluon plasma through the suppression of the back-to-back peak that signals the azimuthal decorrelation of the di-jet system [24,25].

Recently, the possibility of measuring TMB using jet substructure observables like the Soft Drop grooming angle has also been investigated [26,27].

TMB is encoded in the so-called jet quenching parameter $\hat{q}$, which is roughly speaking the average transverse momentum squared $k_{\perp}^2$ acquired per unit of time by a fast parton propagating in a dense medium

$$ \hat{q} \sim \frac{d\langle k_{\perp}^2 \rangle}{dt}. \qquad (1) $$

In the absence of quantum corrections, this simple Brownian diffusion in transverse space picture leads to a TMB distribution which is peaked around a characteristic transverse momentum scale, the saturation scale $Q_s$ such that

$$ Q_s^2 \sim \hat{q}L, \qquad (2) $$

for a given system size $L$. On the other hand, at large $k_{\perp}$, the distribution becomes dominated by rare, single hard scattering with a medium quasiparticle, and therefore displays the characteristic Rutherford-like power law $1/k_{\perp}^4$. The latter has received a lot of attention in the past few years, as it would signal the presence of pointlike quasiparticles in the quark-gluon plasma [28,29]. Of course, the study of this regime in the context of jets in heavy ion collisions is rather challenging due to the large background of soft particles and the interplay with inelastic higher order processes, on the one hand, and to the low cross section associated with these rare events on the other hand.

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.

<footnote>*caucal@subatech.in2p3.fr</footnote>

<footnote>†mehtartani@bnl.gov</footnote>
On the theory side, a lot of progress has been made in the description of TMB. A compact formula à la Molière [30] encompassing both the multiple soft scattering and single hard scattering regimes has been derived in [31]. Furthermore, the effects of transverse flow on the TMB distribution have been investigated analytically and numerically in [32–34]. Also, in the small transverse momentum domain of the distribution, where nonperturbative physics dominates, recent lattice calculations opened up the possibility of achieving a comprehensive picture of TMB at all scales [35,36].

Quantum corrections to TMB have also triggered many recent studies since the seminal papers [37,38] in which the authors show that radiative corrections are enhanced by a double logarithm of the system size L, $\hat{q}_{NLO} \sim \hat{q}_{LO} \times \alpha_s \ln^2 L$. These quantum corrections of order $\alpha_s \propto g^2$ are of a different origin than the classical corrections of order g associated with soft thermal modes in the plasma [39] (for a recent discussion about the interplay between quantum and classical corrections to $\hat{q}$, see [40]). The logarithmic dependence on the system size encoded in higher order corrections is an expression of the nonlocal nature of the quantum fluctuations that were shown to qualitatively affect the underlying diffusive process.

Remarkably, these quantum corrections were shown to be process independent up to single logarithmic accuracy (at large $N_c$), since they also appear in NLO corrections to the medium-induced gluon spectrum [41–44]. This suggests a renormalization group approach in order to resum these potentially large logarithmic terms. Upon resummation to all orders, the resulting TMB distribution showcases interesting physical properties: the typical width of the distribution, parametrically given by the saturation momentum $Q_s$, grows with the system size L faster than the standard diffusive exponent $L^{1/2}$, and the large-$k_\perp$ tail is a power law which deviates from the Rutherford behavior by a $\sqrt{\alpha_s}$ correction [45]. These features are characteristic of anomalous (super) diffusive processes. Moreover, TMB distribution and $Q_s$ tend to universal limits at large system sizes, such that these quantities are not anymore sensitive to the details of the initial conditions and the non-perturbative regime.

So far, the resummation of higher orders of $\hat{q}$ and the TMB distribution have been obtained at double logarithmic accuracy (DLA), meaning that only terms of the form $\alpha_s^n \ln^{2n} L$ are included to all orders in perturbation theory.

In this paper, we address for the first time the effect of single logarithmic corrections on the saturation momentum $Q_s(L)$ in the asymptotic regime where L is large. Relying on the mapping between the evolution equation for $\hat{q}$ and the equations governing the propagation of traveling wave fronts in nonlinear physics [46–50], we can compute the large L asymptotic behavior of the saturation momentum, even without having analytic control of the full nonlinear evolution equation [45,51]. We argue that the evolution of $Q_s$ is dominated by the double logarithmic regime of QCD,

namely by soft and collinear gluon radiations. Therefore, the corrections beyond the DLA can be obtained either from a BFKL [52,53] or a DGLAP [54–56] approach, as the double logarithmic regime is common to these two equations. We choose to proceed using the BFKL equation at leading and next-to-leading logarithmic accuracy (the DGLAP case is presented in Appendix A), as it is more easily justified from a physical point of view and simpler to implement in practice.

The main results of the paper are the universal asymptotic behavior of $Q_s$ at 3-loops in planar conformal $\mathcal{N}=4$ SYM theory and to all loop orders in QCD with massless quarks and at large $N_c$. For QCD, it reads

$$ \ln Q_s^2(L) = Y + 2\sqrt{4b_0Y} + 3\xi_1(4b_0Y)^{1/6} - \frac{1}{4}(3 + 8b_0)\ln Y + \mathcal{O}(1), \quad (3) $$

where $Y = \ln(L)$ and $b_0 = 1/\beta_0$ is the inverse of the 1-loop coefficient of the QCD $\beta$-function. The subleading terms in the expansion (3) can be found in Eq. (58). We refer to this result as an “all-order” result (formally achieved in the asymptotically large $Y = \ln(L)$ limit) as each coefficient in this asymptotic series is exact, and therefore valid to all orders in perturbative QCD. It is then quite remarkable that an exact, all order result can be obtained from perturbative QCD. In $\mathcal{N}=4$ SYM theory, the convergence of the perturbative series seems to be very fast, as we notice a very mild modification from the 2-loop to the 3-loop result. We observe the same feature in QCD: the corrections to the universal asymptotic behavior beyond the DLA results [51,57] are small, and therefore, the DLA with running coupling turns out to be a very good approximation even at moderate values of L.

This paper is organized as follows: in the next section, we briefly review the calculation of the saturation momentum at tree level and 1-loop. We also set up our notations for the rest of the paper. Section III discusses the nonlinear quantum evolution of the dipole cross section and the connection with the BFKL equation in the dilute, linear regime which drives the universality aspects of the saturation momentum for large system sizes. In Sec. IV, we revisit the fixed coupling evolution of the jet quenching parameter from the BFKL language, and address the corrections beyond the double logarithmic approximation in fixed coupling QCD and planar $\mathcal{N}=4$ SYM theory. Finally, in Sec. V, we solve the running coupling evolution at single logarithmic accuracy and present our all order result for the universal behavior of $Q_s(L)$.

## II. THE SATURATION MOMENTUM AT TREE LEVEL AND 1-LOOP

The central object of this paper is the transverse momentum (TMB) distribution $P(k_\perp, t)$ of a high energy parton in the color representation $R=A, F$. It represents
the probability to acquire a transverse momentum $k_{\perp}$ after a time $t$. We focus on the regime in which the momentum transfer from the medium is much smaller than the incoming parton energy $E = P^{+}$ such that the medium interactions do not alter significantly its direction of propagation. Under this approximation, the TMB distribution can be related to the forward scattering amplitude $S(r_{\perp})$ of an effective dipole with transverse size $r_{\perp}$ via a Fourier transform:

$$ \mathcal{P}(k_{\perp}) = \int d^2x_{\perp} e^{-ik_{\perp} \cdot r_{\perp}} S(r_{\perp}). \quad (4) $$

Assuming the interactions between the dipole and the medium scattering centers are local and instantaneous, with a collision rate $C(q_{\perp})$, the forward scattering amplitude exponentiates as follows:

$$ S(r_{\perp}) = \exp \left[ -\frac{C_R}{N_c} \sigma_{\text{dip}}(r_{\perp}) L \right], \quad (5) $$

where $\sigma_{\text{dip}}$ is the so-called dipole cross section. At leading order, it is related to the collision rate $C(q_{\perp})$ according to

$$ \sigma_{\text{dip}}(r_{\perp}) = \int \frac{d^2 q_{\perp}}{(2\pi)^2} (1 - e^{iq_{\perp} \cdot r_{\perp}}) C(q_{\perp}), \quad (6) $$

$$ = \frac{1}{4} \hat{q}(1/r_{\perp}^2, L) r_{\perp}^2 (1 + \mathcal{O}(r_{\perp}^2 \mu^2)). \quad (7) $$

In the second equality, $\mu$ is a nonperturbative infrared momentum scale, typically of order of the plasma Debye mass $m_D \sim gT$ for a plasma at temperature $T$. Equation (7) essentially defines the jet quenching parameter $\hat{q}(k_{\perp}^2, L)$ (in the adjoint representation) in the perturbative regime. At tree level and for a static medium, $\hat{q}$ does not depend on the system size $L$. Applying the gradient $\nabla_{r_{\perp}}$ twice on Eq. (6), and assuming that $\hat{q}$ is a weak function of $r_{\perp}$, it is straightforward to see that it can equivalently be defined as the second moment of the collision rate $C(q_{\perp})$, with an UV momentum cutoff set by $1/r_T$.

The saturation momentum is an emergent scale resulting from the unitarization of the TMB distribution at small $k_T$ (but still much larger than $\mu$). It controls the transition between the dilute regime in which the TMB distribution has the typical Rutherford power law decay $\sim 1/k_{\perp}^4$ and the dense regime where the physics of multiple soft scatterings dominates, typically about $Q_s \gg \mu$. This transition scale is defined by the implicit relation [31,58,59]

$$ S(1/Q_s^2(L)) \equiv e^{-1/4} \Leftrightarrow \hat{q}(Q_s^2(L), L)L \equiv Q_s^2(L). \quad (8) $$

The number $e^{-1/4}$ is arbitrary here, and we will address the sensitivity of our results to this choice later in this paper. At tree level, using the hard thermal loop (HTL) result for the

collision rate $C(q_{\perp})$ that correctly accounts for modes with $|q_{\perp}| \ll T$, one finds that

$$ \hat{q}^{(0)}(k_{\perp}^2, L) = \hat{q}_0 \ln(k_{\perp}^2/\mu^2), \quad (9) $$

with $\hat{q}_0 = \alpha_s N_c m_D^2 T$ and $\mu = m_D e^{-1+\gamma_E}/2$. The exact fixed coupling values for $\hat{q}_0$ and $\mu$ that also includes modes with $|q_{\perp}| \gtrsim T$ can be found in [39,60]. When using the 1-loop running coupling the logarithmic dependence upon the hard scale $|k_{\perp}| \sim 1/|r_{\perp}|$ disappears resulting in $\hat{q}^{(0)}$ being a constant coefficient with $\hat{q}^{(0)} \propto \alpha_s(m_D^2)n$ [61–63], where $n$ is the density of scattering centers. All these details will not matter in the following discussion, owing to the universal property of the asymptotic regime of $Q_s$. Also, note that $\hat{q}_0$ is proportional to $\alpha_s n$. Therefore, this “tree-level” computation is actually an all order resummation in the number of interactions with medium scattering centers. This should be kept in mind since when we will consider the weak coupling limit $\alpha_s \to 0$, it will be implicitly assumed that this limit is taken with the product $\alpha_s nL$ [that appears in the exponential Eq. (5)] fixed.

Using the expression (9), we get the following result for the saturation momentum as a function of $L$:

$$ Q_s^{2,(0)}(L) = \hat{q}_0 L W_{-1} \left( -\frac{\mu^2}{\hat{q}_0 L} \right) \approx \hat{q}_0 L \ln \left( \frac{\hat{q}_0 L}{\mu^2} \right), \quad (10) $$

where $W_p(x)$ is the Lambert function on the $p$th branch. In what follows, it will be convenient to proceed with the following variables:

$$ \rho_s(Y) = \ln(Q_s^2(L)/\mu^2), \quad Y = \ln(L/\tau_0), \quad (11) $$

with $\tau_0 \equiv \mu^2/\hat{q}_0$, the asymptotic behavior of the saturation momentum at tree level is given by

$$ \rho_s^{(0)}(Y) = Y + \ln(Y) + \dots \quad (12) $$

At NLO, the TMB distribution, and therefore the saturation momentum itself are enhanced by large double logarithms $\alpha_s \ln^2(L/\tau_0)$ of the system size. References [37,43] report the following 1-loop result:

$$ Q_s^2(L) = Q_s^{2,(0)}(L) \times \left[ 1 + \bar{\alpha}_s \ln^2 \left( \frac{L}{\tau_0} \right) + \bar{\alpha}_s \left( 2 \ln(2) - \gamma_E - \frac{1}{3} \right) \ln \left( \frac{L}{\tau_0} \right) + \mathcal{O}(\alpha_s) \right], \quad (13) $$

with $\bar{\alpha}_s = \alpha_s N_c / \pi$ and where the $\mathcal{O}(\alpha_s)$ denotes the $\alpha_s$ finite terms. $\gamma_E$ is the Euler-Mascheroni constant. The potentially large double logarithm in Eq. (13) is the dominant radiative correction in the regime $E \gg \omega_c = \hat{q}L^2/2$ we are working, and needs to be resummed to all orders in perturbation theory when $\alpha_s Y^2 = \mathcal{O}(1)$.
### III. NONLINEAR EVOLUTION OF THE DIPOLE CROSS SECTION

A nonlinear evolution equation resumming the double and single logarithms of Eq. (13) has been proposed in [42]. It is formulated directly in terms of the dipole cross section $\sigma_{\text{dip}}(\mathbf{r}_{\perp}, \omega)$ which acquires a rapidity $\ln \omega$ dependence through the evolution, where $\omega \equiv k^{+}$ the light-cone energy of the gluon fluctuation. Schematically, this evolution equation reads

$$ \frac{\partial \sigma_{\text{dip}}}{\partial \ln \omega} = \mathcal{H}_{\hat{q}} \otimes \sigma_{\text{dip}}(\mathbf{r}_{\perp}, \omega), \quad (14) $$

where $\mathcal{H}_{\hat{q}}$ is a nonlinear operator whose precise definition is not important to us [see Eq. (4.24) in [42]]. It satisfies the property that after one step of the evolution Eq. (14), one gets the double and single logarithmic terms in Eq. (13).

Such evolution equation is difficult to solve both analytically and numerically. However, in this paper, we are mainly interested in the large L limit of the saturation momentum $Q_s(L)$ arising from Eq. (14). In a series of recent papers, a new mathematical method based on the analogy between the evolution equation (14) and the propagation of traveling wave fronts into unstable states has been developed in order to compute analytically the asymptotic behavior of $Q_s(L)$, despite the absence of general analytic solutions of Eq. (14). In particular, the existence of traveling wave solutions to Eq. (14) allows us to simplify the problem and study the dilute (or linear) regime which drives the growth of the perturbations around the unstable state that determines the speed of the front $\dot{\rho}_s = d\rho_s/dY$, in the presence of an absorptive boundary at $Q_s$. Essentially, the details of nonlinear dynamics responsible for the saturation (unitarization) of the dipole S matrix at $Q_s$ are irrelevant for the determination of the universal behavior of the saturation scale [50,64] (see also [48,65–68] for concrete examples in the context of small x physics).

In order to define the dilute regime, we first recall that the formation time of the gluon fluctuation is of order $\tau \equiv 1/k^{-} = 2\omega/k_{g\perp}^2$. In the dilute regime, the transverse momentum accumulated via multiple collisions over the formation time of the gluon, of order $\hat{q}\tau$, cannot exceed its transverse momentum. This criterion guarantees that the gluon is not sensitive to multiple soft scatterings over its formation time, so that only a single (or few [40], near the boundary of the constraint) scattering contributes to the cross section. Therefore, the dilute regime corresponds to

$$ k_{g\perp}^2 \gg Q_s^2(\tau) \sim \hat{q}\tau, \quad (15) $$

or, in terms of $\omega$, $k_{g\perp}^2 \gg \sqrt{\hat{q}\omega}$. In this regime, the nonlinear evolution equation reduces to the well-known BFKL equation [42]

$$ \frac{\partial \sigma_{\text{dip}}(\mathbf{r}_{\perp}, \omega)}{\partial \ln \omega} = \frac{\bar{\alpha}_s}{2\pi} \int d^2\mathbf{z}_{\perp} \mathcal{K}_{r0z}[\sigma_{\text{dip}}(\mathbf{r}_{\perp} - \mathbf{z}_{\perp}, \omega) + \sigma_{\text{dip}}(\mathbf{z}_{\perp}, \omega) - \sigma_{\text{dip}}(\mathbf{r}_{\perp}, \omega)], \quad (16) $$

with the LO BFKL kernel in coordinate space

$$ \mathcal{K}_{xyz} = \frac{(\mathbf{x}_{\perp} - \mathbf{y}_{\perp})^2}{(\mathbf{z}_{\perp} - \mathbf{x}_{\perp})^2 (\mathbf{z}_{\perp} - \mathbf{y}_{\perp})^2}, \quad (17) $$

and $\bar{\alpha}_s \equiv \alpha_s C_A/\pi$, albeit with an additional constraint that enforces the condition Eq. (15) in coordinate space [see for instance Eq. (23) below]. We shall discuss the scale choice for the running coupling in the case of the jet quenching problem in the next section. In the BFKL equation in coordinate space (16), we have parametrically $|\mathbf{k}_{\perp}| \sim 1/|\mathbf{r}_{\perp}|$, (we remind that $\mathbf{k}_{\perp}$ is the final transverse momentum of the incoming parton) and $|\mathbf{k}_{g\perp}| \sim 1/|\mathbf{z}_{\perp}|$ for the transverse momenta of the gluons along the ladder.

In addition to [42], the relevance of the BFKL equation in the renormalization of the jet quenching parameter has also been discussed from an effective field theory perspective in [69,70]. However, as we have anticipated it, the dynamics being of double logarithmic nature BFKL and DGLAP evolutions are equally good to compute the asymptotics of the saturation scale for jet quenching. We stress that this is quite different in studies of proton structure at small x where the rapidity logarithm, $Y \equiv \ln 1/x$, is assumed to be much larger than the collinear logarithm $\ln k_{\perp}^2 < \ln Q^2 \ll \ln 1/x$. We will return to this crucial difference between small-x and jet quenching evolution, and its implications, when we discuss the BFKL kernel in Mellin space.

If one aims at reaching single logarithmic accuracy for $\sigma_{\text{dip}}(\mathbf{r}_{\perp}, \omega)$, one needs also to include the NLL BFKL evolution. This contribution is not manifest in the fixed order computation in [37] because it appears at 2-loops in perturbation theory (it is a corrections of order $\alpha_s^2$). However, it is accompanied by a double log $\sim Y^2$, and therefore matters at single log accuracy since $\alpha_s^2 Y^2 = \mathcal{O}(1)$.

It is convenient to write the full NLL BKFL equation using the variable $\rho = -\ln(\mathbf{r}_{\perp}^2\mu^2)$ and $\eta = \ln(\omega/\omega_0)$ with the infrared energy scale $\omega_0 = \hat{q}_0\tau_0^2/2$,

$$ \frac{\partial \sigma_{\text{dip}}}{\partial \eta} = [\bar{\alpha}_s \chi_{\text{LL}}(-\partial_{\rho}) + \bar{\alpha}_s^2 \chi_{\text{NLL}}(-\partial_{\rho})] \sigma_{\text{dip}}(\rho, \eta). \quad (18) $$

In this equation, $\chi_{\text{LL}}$ and $\chi_{\text{NLL}}$ are the Mellin space representation of the BFKL leading log and next-to-leading log kernel. $\chi_{\text{LL}}$ has the familiar expression

$$ \chi_{\text{LL}}(\gamma) = 2\psi(1) - \psi(\gamma) - \psi(1-\gamma), \quad (19) $$

with $\psi(x)$ the digamma function. The expression for $\chi_{\text{NLL}}$ is more complicated and can be found in [71]. In what
follows, only the pole structure of $\chi_{\text{NLL}}$ at $\gamma = 1$ matters. Indeed, the jet quenching evolution problem is intrinsically double logarithmic in nature, in the sense that it is dominated by gluon fluctuations which are both strongly ordered in transverse momenta in the collinear regime $k_{\perp}^2 \gg k_{g\perp} \gg \dots \gg \mu^2$ and strongly ordered in energy or light-cone plus component $E \gg \omega \gg \dots \gg \omega_0$. This double logarithmic regime is driven by the poles in $\gamma = 1$ of $\chi_{\text{LL}}$ and $\chi_{\text{NLL}}$. In contrast, the standard small-$x$ evolution is driven by finite values of $\gamma$, i.e., $1/2$ and $0.372$ for BFKL with and without a saturation boundary, respectively [48,65,72,73].

Before proceeding further we need first to address the issue of the large NLL correction in BFKL equation that are associated with a spurious triple pole in $\chi_{\text{NLL}}$ at $\gamma = 1$ due to the wrong choice of the evolution variable $\omega$ in the collinear regime [74–77]. The leading behavior of $\chi_{\text{NLL}}$ is indeed [71]

$$ \chi_{\text{NLL}}(\gamma) = -\frac{1}{(1-\gamma)^3} + \frac{B_g}{(1-\gamma)^2} + \mathcal{O}\left(\frac{1}{1-\gamma}\right), \quad (20) $$

with

$$ B_g = -\frac{11}{12} - \frac{N_f}{6N_c^3} \approx -\frac{11}{12}. \quad (21) $$

In the second equality, we have used the large $N_c$ approximation. The solution to this issue is well known. Instead of using $\omega = k^+$ as the evolution variable, one should use the lifetime $\tau = 1/k^- = 2\omega/k_{g\perp}^2$ of the gluon fluctuation. In terms of $Y = \ln(\tau/\tau_0)$, the NLL evolution equation for $\sigma_{\text{dip}}(\rho, Y)$ is identical to Eq. (18) (up to pure $\alpha_s$ corrections), but this time, $\chi_{\text{NLL}}$ has no triple poles anymore [77,78]. The equation we shall study is therefore

$$ \frac{\partial \sigma_{\text{dip}}}{\partial Y} = [\bar{\alpha}_s \chi_{\text{LL}}(-\partial_\rho) + \bar{\alpha}_s^2 \tilde{\chi}_{\text{NLL}}(-\partial_\rho)] \sigma_{\text{dip}}(\rho, Y), \quad (22) $$

with $\tilde{\chi}_{\text{NLL}} = \chi_{\text{NLL}} + 1/(1-\gamma)^3 + \mathcal{O}(1)$. Note that this equation has other issues in the double logarithmic anticolinear regime [78], but this regime is irrelevant in our case, as we shall see in the next section.

Finally, even though the dilute (linear) regime of the nonlinear evolution drives the asymptotic behavior of $\rho_s$, it is important to keep in mind that there is a major difference between BFKL evolution and the problem at hand. The latter pertains to the existence of a saturation (absorptive) boundary in the double logarithmic phase space. This saturation boundary can be accounted for via the following step function in coordinate space representation of the BFKL equation [using $1/k_{g\perp}^2 \sim \max((x_\perp - z_\perp)^2, z_\perp^2)$]:

$$ \Theta\left(\frac{4}{Q_s^2(\tau)} - \max((x_\perp - z_\perp)^2, z_\perp^2)\right), \quad (23) $$

where $Q_s^2(\tau) = \hat{q}\tau$ is the saturation momentum evaluated at the gluon formation time $\tau \equiv 1/k^-$. This step function enforces the transverse momentum of the gluon fluctuation to be larger than the saturation momentum so as the fluctuation is not affected by the LPM effect.

It is worth noting that because the dynamics is dominated by strongly ordered transverse sizes in DLA, i.e., $x_\perp \ll z_\perp$ the theta function can be simplified as $\Theta(4/Q_s^2(\tau) - z_\perp^2)$. With this simplification, one easily checks that Eq. (16) with this constraint reduces to the nonlinear evolution equation for $\hat{q}$ written in [41,42,45] [see also Eq. (42) below].

We are now left with the study of the NLL BFKL equation with a saturation boundary (or LO DGLAP as discussed in Appendix A). The exact implementation of this saturation condition is not decisive in the asymptotic regime, however, its very existence constrains the shape of the traveling wave ansatz that we shall discuss in what follows.

## IV. REVISITING THE FIXED COUPLING PROBLEM

In order to make the connection with our previous studies in [45,51], we shall consider the evolution of $\hat{q}(\rho, Y)$ instead of $\sigma_{\text{dip}}(\rho, Y)$. The function $\hat{q}(\rho, Y)$ is defined according to Eq. (7), or equivalently in terms of $\rho$,

$$ \hat{q}(\rho, Y) = 4\mu^2 e^\rho \sigma_{\text{dip}}(\rho, Y). \quad (24) $$

Plugging this definition into Eq. (22), using $e^\rho\chi(-\partial_\rho)e^{-\rho} = \chi(1-\partial_\rho)$ and the symmetry property $\chi(\gamma) = \chi(1-\gamma)$ of the kernel, one ends up with the following equation for $\hat{q}$

$$ \frac{\partial \hat{q}(\rho, Y)}{\partial Y} = [\bar{\alpha}_s \chi_{\text{LL}}(\partial_\rho) + \bar{\alpha}_s^2 \tilde{\chi}_{\text{NLL}}(\partial_\rho)] \hat{q}(\rho, Y). \quad (25) $$

The saturation condition Eq. (23), that translates into $\rho > \rho_s(Y)$ in the new variables, is implicit in the above equation.

The purpose of this section is to first recover the known results [51,57] for the asymptotic expansion of $\rho_s$ at DLA and fixed coupling from this equation and to discuss the case of the conformal $N=4$ SYM theory in which the coupling does not run. We also demonstrate that contrary to the BFKL problem which is driven by the behavior of the kernel around $\gamma_c \approx 0.327$, the evolution of $\hat{q}$ is governed by the double logarithmic collinear regime and therefore driven by the behavior of the kernel around $\gamma = 0$ (or $\gamma = 1$ for the dipole cross section) as stated in the previous section.

At DLA and fixed coupling, one can neglect the NLL term and simply have

$$ \frac{\partial \hat{q}(\rho, Y)}{\partial Y} = \bar{\alpha}_s \chi_{\text{LL}}(\partial_\rho) \hat{q}(\rho, Y), \quad (26) $$
for a fixed $\bar{\alpha}_s$. The key starting point of our analysis consists of using the ansatz

$$ \hat{q}(\rho, Y) = e^{\rho_s(Y)-Y} e^{\beta x} f(x, Y), \quad x = \rho - \rho_s(Y) \quad (27) $$

in order to solve Eq. (26) perturbatively in the limit $Y \to \infty$. It is motivated by the presence of the nonlinear saturation condition in the evolution, $Q_s^2(L) = \hat{q}(Q_s^2(L), L)L$ which reads

$$ \hat{q}(\rho_s(Y), Y) = e^{\rho_s(Y)-Y} \quad (28) $$

in terms of $\rho_s$ and $Y$. Therefore, the function $f$ satisfies $f(0, Y) = 1$ for all $Y$. This identity (more precisely, the 1 on the r.h.s.) depends on the definition we adopt for $Q_s$.

The form of this ansatz, though inspired by the study of traveling waves propagation governed by the Balitsky-Kovchekov (BK) equation [79,80] in small $x$ physics, differs from the latter in a crucial way. It involves an additional kinematic factor $1/\tau \sim e^{-Y}$ which is responsible for driving the evolution towards the double logarithmic regime $\gamma \sim 0 + \mathcal{O}(\bar{\alpha}_s^{1/2})$ instead of $\gamma_c \sim 0.327$. Nevertheless, we will see that our BFKL equation with saturation constraint admits traveling wave solutions (or geometric scaling solutions) of the form $\hat{q}(Y, \rho) = e^{\rho_s(Y)-Y} e^{\beta x} f(x)$ in the limit $Y \to \infty$. Since $x = \rho - \rho_s(Y)$, such solutions correspond indeed to the propagation of a front at the speed $d\rho_s/dY$ from the left to the right along the $\rho$ axis.

As alluded to below Eq. (8), the saturation momentum is to be defined up to a undetermined multiplicative constant. This freedom can be absorbed into a redefinition of $Y$ through a constant shift (equivalent to a redefinition of the nonperturbative scale $\tau_0$). Such redefinition does not affect the universal terms in the asymptotic development of $\rho_s$ that we intend to calculate here, as can be checked by shifting $Y$ in the expressions (46) and (58).

From the existence of a scaling limit [45], we expect the function $f(x, Y)$ to converge towards a well-defined function $f(x)$ as $Y \to \infty$. The parameter $\beta$ will be determined later. Plugging this ansatz into Eq. (26), and expanding the kernel $\chi_{LL}(\gamma)$ around $\gamma = \beta$, we find

$$ (\dot{\rho}_s - 1 - \dot{\rho}_s \beta) f - \dot{\rho}_s \partial_x f + \partial_Y f = \bar{\alpha}_s \sum_{p=0}^{\infty} \frac{\chi_{LL}^{(p)}(\beta)}{p!} \partial_x^p f, \quad (29) $$

where

$$ \dot{\rho}_s(Y) \equiv \frac{d\rho_s(Y)}{dY} \quad \text{and} \quad \chi_{LL}^{(p)}(\beta) \equiv \frac{\partial^p \chi_{LL}(\beta)}{\partial \beta^p}. \quad (30) $$

Taking the limit $Y \to \infty$ and using the existence of a scaling limit for $f$, we end up with the following relations between the speed of the front

$$ c = \lim_{Y \to \infty} \dot{\rho}_s(Y), \quad (31) $$

and the critical value $\beta_c$ of $\beta$ that minimizes [45] the velocity $c$,

$$ c - 1 - c\beta_c = \bar{\alpha}_s\chi_{LL}(\beta_c), \quad (32) $$

$$ -c = \bar{\alpha}_s \chi'_{LL}(\beta_c), \quad (33) $$

by simply identifying the terms proportional to $f$ and $\partial_x f$. This system cannot be solved analytically, however one can find the series expansion of $c$ and $\beta_c$ in powers of $\bar{\alpha}_s$:

$$ c = 1 + 2\sqrt{\bar{\alpha}_s} + 2\bar{\alpha}_s + \mathcal{O}(\bar{\alpha}_s^{3/2}), \quad (34) $$

$$ \beta_c = \sqrt{\bar{\alpha}_s} - \bar{\alpha}_s + \mathcal{O}(\bar{\alpha}_s^{3/2}). \quad (35) $$

At this stage, only the terms of order up to $\mathcal{O}(\sqrt{\bar{\alpha}_s})$ are under control since we neglect the NLL BFKL term in the evolution. To recover the DLA results reported in [45], one can approximate $\chi_{LL}$ by its most singular behavior as $\gamma \to 0$. The LL BFKL kernel behaves like

$$ \chi_{LL}(\gamma) = \frac{1}{\gamma} + 2\zeta(3)\gamma^2 + \mathcal{O}(\gamma^4), \quad (36) $$

at small $\gamma$, with $\zeta(x)$ the Riemann zeta function and $\zeta(3) \approx 1.2$. Using $\chi_{LL} = 1/\gamma$, one can solve exactly the system (34)–(35) (assuming $c \ge 1$) and find $c = 1 + 2\sqrt{\bar{\alpha}_s + \bar{\alpha}_s^2} + 2\bar{\alpha}_s$ and $\beta_c = \sqrt{\bar{\alpha}_s + \bar{\alpha}_s^2} - \bar{\alpha}_s$ in agreement with [45]. The $\mathcal{O}(\sqrt{\bar{\alpha}_s})$ and $\mathcal{O}(\alpha_s)$ terms are not affected by the use of the full BFKL kernel or its approximation $\chi_{LL} \sim 1/\gamma$ (contrary to the BK case in small $x$ physics). One can actually check that the deviations enter at $\mathcal{O}(\alpha_s^2)$:

$$ \beta_c = \bar{\alpha}_s^{1/2} - \bar{\alpha}_s + \frac{1}{2}\bar{\alpha}_s^{3/2} + 2\zeta(3)\bar{\alpha}_s^2 + \mathcal{O}(\bar{\alpha}_s^{5/2}). \quad (37) $$

The last term cannot be obtained from the approximated kernel, as can be checked from the expansion of $\sqrt{\bar{\alpha}_s + \bar{\alpha}_s^2} - \bar{\alpha}_s$. It is a contribution beyond the double logarithmic regime of BFKL which would enter at the same order as the most singular term of the N$^3$LO BFKL kernel. This demonstrates that the asymptotic behavior of $\dot{\rho}_s$ and $\hat{q}(\rho, Y)$ is mainly sensitive to the collinear double logarithmic regime of the BFKL evolution, since the details of the leading log kernel $\chi_{LL}$ are not important up to order $\alpha_s^2$ (4-loops) in pQCD. In other words, the saddle point $\beta_c \sim \sqrt{\bar{\alpha}_s}$ is close to 0 for the $\hat{q}$ evolution, contrary to the case of small-$x$ evolution of the dipole operator where the saddle point lies at $\beta_c \approx 0.6275$.

Since the physics is dominated by the double logarithmic regime, which is common to both BFKL and DGLAP, one should be able to recover the results obtained in this paper
<figure><img src="image_3.png" /><figcaption>FIG. 1. The plane $(\rho, Y)$ and the standard DGLAP and BFKL regime of pQCD. The saturation scale $Q_s$ evolves along the diagonal of this diagram since $\rho \sim Y$.</figcaption></figure>

from a BFKL approach using instead a DGLAP-like evolution with $\rho$ as the evolution variable. It turns out to be the case as shown in Appendix A. This is also illustrated in Fig. 1 where we display the $(\rho, Y)$ domain of pQCD as well as the standard DGLAP and BFKL directions along which the gluon distribution is evolved. The evolution of $Q_s$ follows the diagonal of this diagram as a consequence of the constraint $\rho \sim \rho_s \sim Y$, which is also where the DGLAP and BFKL evolution “merge” in the double logarithmic approximation. It is therefore natural that the corrections beyond DLA can be obtained from both DGLAP or BFKL approaches.

If one includes the NLL BFKL kernel $\tilde{\chi}_{\text{NLL}}$, the system of Eqs. (34)–(35) is modified and the right-hand side receives a contribution $\tilde{\alpha}_s^2 \tilde{\chi}_{\text{NLL}}(\beta_c)$:

$$c - 1 - c\beta_c = \bar{\alpha}_s\chi_{\text{LL}}(\beta_c) + \bar{\alpha}_s^2\tilde{\chi}_{\text{NLL}}(\beta_c), \quad (38)$$

$$-c = \bar{\alpha}_s\chi'_{\text{LL}}(\beta_c) + \bar{\alpha}_s^2\tilde{\chi}'_{\text{NLL}}(\beta_c). \quad (39)$$

With this additional term, we gain control over the terms of order $O(\bar{\alpha}_s)$ in the $\alpha_s$ expansion of $c$ and $\beta_c$. Again, since the problem is dominated by the singular behavior at $\gamma = 0$, it is sufficient to use the approximation $\tilde{\chi}_{\text{NLL}}(\gamma) = B_g/\gamma^2$, and one finds that

$$c = 1 + 2\sqrt{\bar{\alpha}_s} + (2 + B_g)\bar{\alpha}_s + O(\bar{\alpha}_s^{3/2}). \quad (40)$$

This result extends to single logarithmic accuracy the value of the traveling wave speed $c$.

The term of order $O(\bar{\alpha}_s^{3/2})$ receives contribution from both the subleading pole in $1/\gamma$ in the NLL BFKL kernel and from the pole in $1/\gamma^3$ of the N$^2$LL BFKL kernel. In fact, since the pole structure of the NLL and N$^2$LL BFKL equations are known thanks to the DGLAP/BFKL

duality [81,82], $\tilde{\chi}_{\text{NLL}} \sim B_g/\gamma^2 + a_{1,-1}/\gamma$ and $\tilde{\chi}_{\text{N$^2$LL}} \sim a_{2,-3}/\gamma^3$, the value of $c$ can be known up to order $\bar{\alpha}_s^{3/2}$, with the coefficient of the $\bar{\alpha}_s^{3/2}$ term equal to $1 + 3B_g - B_g^2 + a_{1,-1} + a_{2,-3}$.

As an illustrative example, let us consider the case of planar $\mathcal{N} = 4$ SYM theory, for which there are many results on the BFKL regime. In particular, we rely on the calculation of the 3-loops N$^2$LL BFKL/BK equation in this theory [83–85]. The “classical” $O(g)$ corrections to $\hat{q}$ are also known in this theory [86], but such corrections should be included in the initial condition of the quantum evolution, with proper matching [40]. In planar $\mathcal{N} = 4$ SYM theory, since the coupling does not run ($B_g = 0$) and $a_{1,-1} = 0$ [87], we have the expansion

$$c = 1 + 2\sqrt{\bar{\alpha}_s} + 2\bar{\alpha}_s + (1 + a_{2,-3}^{\mathcal{N}=4})\bar{\alpha}_s^{3/2} + O(\bar{\alpha}_s^2), \quad (41)$$

with the coefficient of the triple pole of the 3-loops BFKL kernel $a_{2,-3}^{\mathcal{N}=4} = -\zeta(2) = -\pi^2/6$ [83–85]. The behavior of this series is shown Fig. 2, and compared to the strong coupling limit $c = 2$ obtained from AdS/CFT [88–90].

It would be interesting to exploit the BFKL/DGLAP duality (or existing results for the N$^p$LL BFKL equation with $p \ge 4$ [91]) in order to extend this calculation to higher orders in $\alpha_s$ and see the convergence of the series and the approach towards the strong coupling regime.

Note that although this is a rather academic exercise given that the running of the coupling spoils these conclusions, as we shall see in the next section, the above discussion highlights the fact that in a conformal theory with scale invariant coupling, the calculation of the $\alpha_s^{p/2}$ order of the “anomalous dimension” $c$ of the saturation scale for large media requires the knowledge of the N$^p$LL BFKL equation because the $\alpha_s$ expansion in the fixed

<figure><img src="image_4.png" /><figcaption>FIG. 2. The asymptotic limit of the front wave velocity as a function of the coupling constant $\bar{\alpha}_s = \alpha_s N_c / \pi$ in (planar) $\mathcal{N} = 4$ SYM theory.</figcaption></figure>
coupling case is decoupled from the asymptotic Y expansion. This no longer the case with the running coupling where $\alpha_s$ explicitly depends on Y through its dependence on $Q_s$ and thus, the two expansions are related to one another.

Turning back to QCD, one may want to understand physically where the coefficient of the $\alpha_s$ term in c given by Eq. (40) comes from. The $B_g$ term is the finite part of the gluon splitting function, and therefore this contribution is associated with collinear but nonsoft splittings in the evolution. On the other hand, the 2 term seems more mysterious at first glance. We argue that it comes from the feedback of the quantum evolution of $Q_s$ on the evolution of the dipole cross section $\sigma_{\text{dip}}$ (recall that these physical quantities are related to one another by the nonlinearity of the evolution). Indeed, the fixed coupling DLA equation with $\chi_{\text{LL}} = 1/\gamma$ is equivalent to the equation

$$ \frac{\partial \hat{q}}{\partial Y} = \bar{\alpha}_s \int d^\rho \hat{\rho}(\hat{\rho}', Y), \qquad (42) $$

owing to the formal relation $\int d^\rho d\rho' = 1/\partial_\rho$. In this equation, the lower bound is set by the single scattering criterion $k_{g\perp}^2 \ge Q_s^2$ i.e. $\rho' > \rho_s(Y)$. If we neglect quantum evolution in this lower bound and use the tree level “classical” expression $\rho_s(Y) = Y$ instead, to constraint the $\rho'$ integral such that $\rho' > Y$, this equation can be solved exactly. Defining $Q_s^2(L) = \hat{q}(\hat{q}_0 L, L)L$ in agreement with the previous approximation, the resummed value of the saturation momentum is given by [57,92]

$$ Q_s^2(L) = Q_s^{2,(0)}(L) \times \frac{I_1(2\sqrt{\bar{\alpha}_s Y^2})}{\sqrt{\bar{\alpha}_s Y^2}}, \qquad (43) $$

for a constant initial condition (hence $Q_s^{2,(0)}(L) = \hat{q}_0 L$ here). Here $I_n(x)$ is the modified Bessel function of rank n. This result shows very clearly the double logarithmic resummation structure, since the saturation momentum is expressed as a function of $\alpha_s Y^2 = \alpha_s \ln^2(L/\tau_0)$. In the asymptotic limit, assuming the strong condition $\alpha_s Y^2 \gg 1$ (instead of $\alpha_s Y^2 \sim 1$ only), one can easily derive the behavior of the saturation scale from Eq. (43):

$$ \rho_s(Y) = (1 + 2\sqrt{\bar{\alpha}_s})Y - \frac{3}{2}\ln(Y) + \mathcal{O}(1). \qquad (44) $$

In DLA evolution with linearization of the saturation boundary, one observes that the traveling wave speed is $c = 1 + 2\sqrt{\bar{\alpha}_s}$ without terms of order $\alpha_s$. Hence, one can interpret the 2 term in the $\bar{\alpha}_s$ coefficient of Eq. (40) as a single log effect coming from the quantum evolution of $Q_s$ itself. A similar argument is presented in [57].

In Eq. (44), we also show the subasymptotic correction to the constant speed motion of the saturation front in the

case of the linearized DLA evolution. For the full NLL problem, the subasymptotic corrections to $\dot{\rho}_s$ can be obtained using the same method as in [45,51]. Namely, we expand $f(x, Y)$ as an infinite series in powers of $Y^{-1/2}$ multiplied by diffusive scaling functions $G_n$ [49,72,93]:

$$ f(x, Y) = \sum_{n=-1}^{\infty} Y^{-n/2} G_n \left( \frac{x}{Y^{1/2}} \right). \qquad (45) $$

Using this expansion and solving order by order for the functions $G_n$ with appropriate boundary conditions, we can fix the asymptotic development of $\dot{\rho}_s$. This is shown in Appendix B. We find that this development can be expressed in terms of c and $\beta_c$ as

$$ \dot{\rho}_s = c + \frac{3}{2(\beta_c - 1)} \frac{1}{Y} + \frac{3}{2(\beta_c - 1)^2} \sqrt{\frac{2\pi}{\chi''(\beta_c)}} \frac{1}{Y^{3/2}} + \mathcal{O}(Y^{-2}), \qquad (46) $$

with $\chi''(\beta_c) = \bar{\alpha}_s \chi''_{\text{LL}}(\beta_c) + \bar{\alpha}_s^2 \chi''_{\text{NLL}}(\beta_c) + \dots$. This expression encompasses all universal terms in the expansion of $\dot{\rho}_s$ at fixed coupling. The first nonuniversal terms that are sensitive to the initial condition appear at the order $Y^{-2}$. Such terms can be easily produced by shifting the value of Y by a constant to absorb a change in the nonperturbative parameters $\tau_0$ or $\mu$.

When using the DLA, i.e. $\chi = \bar{\alpha}_s/\gamma$, one recovers from Eq. (46) the expression found in [45,51].

Each coefficient of the Y powers has an $\alpha_s$ expansion which can be obtained from the one of c and the shape of the BFKL kernels near the double logarithmic regime $\gamma = 0$. These expansions are given in Appendix C up to

<figure><img src="image_5.png" /><figcaption>FIG. 3. The front wave velocity as a function of Y in (planar) $\mathcal{N} = 4$ SYM theory with $\bar{\alpha}_s = 0.1$. The red band shows the uncertainty related to the unknown $Y^{-2}$ corrections in $\dot{\rho}_s$, obtained by varying $\tau_0$ by a factor of 2 around a central value. For $Y \gtrsim 5$, next-to-double log corrections dominate over the nonuniversal subasymptotic corrections in $Y^{-2}$.</figcaption></figure>
order $\alpha_s^{3/2}$ in $\mathcal{N} = 4$ SYM theory. They are systemically improvable by including higher logarithmic orders in the BFKL equation. As we shall see in the next section, conformal symmetry breaking in QCD changes this simple picture in a quite dramatic way: the coefficients of the universal terms in the $\rho_s$ expansion are entirely fixed by the BFKL kernel at NLL order.

The expansion (46) of $\dot{\rho}_s(Y)$ as a function of $Y$ is shown in Fig. 3 in the case of $\mathcal{N} = 4$ SYM theory, at DLA, order $\alpha_s$ and order $\alpha_s^{3/2}$ for $\bar{\alpha}_s = 0.1$. The convergence is very good, and one notices that the corrections beyond the double logarithms are more important than the nonuniversal subasymptotic $1/Y^2$ correction (estimated by varying $\tau_0$ by factors between 0.5 and 2) for $Y \gtrsim 4$.

## V. RUNNING COUPLING EFFECTS AND ALL-ORDER RESULT FOR THE SATURATION MOMENTUM

We now address the running coupling problem and compute the universal asymptotic expansion of $\dot{\rho}_s$ in QCD at single logarithmic accuracy. To do so, it is sufficient to consider the 1-loop running coupling:

$$ \alpha_s(\rho) = \frac{b_0}{\rho + \rho_0}, \qquad (47) $$

with $b_0 = 1/\beta_0 = -1/B_g$ and $\rho_0 = \ln(\mu^2/\Lambda_{QCD}^2)$. When including running coupling effects, one has to be careful with the prescription for the running scale. At DLA, it has been argued in [37,57] that the running coupling should be implemented in the following way:

$$ \frac{\partial \hat{q}}{\partial Y} = \int^{\rho} d\rho' \bar{\alpha}_s(\rho') \hat{q}(\rho', Y), \qquad (48) $$

so that the standard DGLAP equation in the DLA is recovered (cf. Appendix A). In Mellin space, the corresponding NLL equation is

$$ \frac{\partial \hat{q}}{\partial Y} = \chi_{LL}(\partial_{\rho})[\bar{\alpha}_s(\rho)\hat{q}(\rho, Y)] + \bar{\alpha}_s^2(\rho)\tilde{\chi}_{NLL}(\partial_{\rho})\hat{q}(\rho, Y), \qquad (49) $$

which is slightly different from the running coupling BFKL equation, where $\bar{\alpha}_s(\rho)$ is outside of the BFKL kernel. In fact, these two prescriptions, or scheme choices, are related to one another by a modification of the NLL kernel [81]. For the NLL term proportional to $\bar{\alpha}_s^2$, the differences between these two schemes are of higher order (NNLL), which explains why we have written $\bar{\alpha}_s^2(\rho)$ outside of the NLL kernel.

We therefore need to study the evolution equation (49) with the pole structure of the NLL BFKL kernel given by Eq. (20). As shown in [51], the DLA running coupling involves a modified geometric scaling limit which takes the form

$$ \hat{q}(\rho, Y) = e^{\rho_s(Y)-Y} e^{\beta_s f(x, Y)}, \qquad x = \frac{\rho - \rho_s(Y)}{\sqrt{Y}}, \qquad (50) $$

with the function $f(x, Y)$ having a scaling limit $f(x)$ as $Y \to \infty$. We then plug this ansatz inside the evolution equation and expand the kernels $\chi_{LL}$ and $\chi_{NLL}$ around the DLA “saddle point” $\beta/\sqrt{Y}$, as expected from the fixed coupling result where $\beta_c \sim \sqrt{\alpha_s(\rho_s)} \propto Y^{-1/2}$. After these manipulations, we find the equation

$$ \begin{aligned} & \left(\dot{\rho}_s - 1 - \frac{\beta_x}{2Y} - \frac{\beta \dot{\rho}_s}{Y^{1/2}}\right) f - \left(\frac{x}{2Y} + \frac{\dot{\rho}_s}{Y^{1/2}}\right) \frac{\partial f}{\partial x} + \frac{\partial f}{\partial Y} \\ & = \sum_{p=0}^{\infty} \frac{\chi_{LL}^{(p)}(\beta/Y^{1/2})}{p! Y^{p/2}} \partial_x^p \left[ \frac{b_0}{\rho_s + xY^{1/2} + \rho_0} f \right] \\ & \quad + \left( \frac{b_0}{\rho_s + xY^{1/2} + \rho_0} \right)^2 \sum_{p=0}^{\infty} \frac{\tilde{\chi}_{NLL}^{(p)}(\beta/Y^{1/2})}{p! Y^{p/2}} \partial_x^p f. \end{aligned} \qquad (51) $$

As in the fixed coupling problem, let us first determine the exact location of the saddle point proportional to $\beta$ and the value of the first nontrivial correction to $\dot{\rho}_s = 1$ such that

$$ \dot{\rho}_s = 1 + \frac{c}{Y^{1/2}} + \dots \qquad (52) $$

This behavior is dictated by the homogeneity of Eq. (51). Seeking for the leading power in $1/\sqrt{Y}$, one gets the two following equations:

$$ c - \beta_c = b_0/\beta_c, \qquad (53) $$

$$ -1 = -b_0/\beta_c^2, \qquad (54) $$

after identification of the coefficients in front of $f$ and $\partial_x f$. Therefore, one recovers the scaling limit of the running coupling evolution equation at DLA:

$$ c = 2\sqrt{b_0}, \qquad (55) $$

$$ \beta_c = \sqrt{b_0}. \qquad (56) $$

We emphasize that these relations are exact to all orders in pQCD. Indeed, contrary to the fixed coupling problem, the NLL kernel (and the higher orders) do not affect the value of $c$ and $\beta_c$. This is a consequence of (i) considering the large system size limit, in which the typical transverse momentum scale in the running coupling is $Q_s \sim \sqrt{\hat{q}_0 L}$ itself and (ii) the asymptotic freedom property of QCD which imposes $\alpha_s$ to decay as $\ln^{-1}(Q_s^2/\Lambda_{QCD}^2) \sim 1/Y$ at large $Y$. QCD conformal symmetry breaking through the running of $\alpha_s$ affects the powers of $Y$ appearing in the asymptotic expansion of $\dot{\rho}_s(Y)$ [and $\rho_s(Y)$ after integration], but the coefficient of the leading power is fully set by the LL BFKL equation, in sharp contrast with the conformal planar $\mathcal{N} = 4$ SYM case.
To obtain the subasymptotic corrections, we use again the leading edge expansion

$$f(x, Y) = \sum_{n=-1}^{\infty} Y^{-n/6} G_n \left( \frac{x}{Y^{1/6}} \right), \quad (57)$$

with a 1/6 diffusive exponent [51]. The interplay between this expansion and the asymptotic series of $\dot{\rho}_s$ enables us to compute the subasymptotic corrections. Namely, the determination of the function $G_{-1}$ fixes the correction of order $Y^{-5/6}$ in $\dot{\rho}_s$, the function $G_0$ determines the correction of order $Y^{-1}$, and so forth (see Appendix B). After a straightforward calculation of these functions similar to the one in [51], we find for the universal asymptotic series,

$$\begin{aligned} \dot{\rho}_s = & 1 + \frac{4b_0}{\bar{Y}^{1/2}} + \frac{2\xi_1 b_0}{\bar{Y}^{5/6}} + b_0(1 - 8b_0 + 4b_0 B_g) \frac{1}{\bar{Y}} \\ & - \frac{7\xi_1^2 b_0}{270 \bar{Y}^{7/6}} - (5 + 1944b_0) \frac{\xi_1 b_0}{81 \bar{Y}^{4/3}} \\ & - 2b_0^2(1 - 8b_0 + 4b_0 B_g) \frac{\ln(\bar{Y})}{\bar{Y}^{3/2}} + \mathcal{O}(\bar{Y}^{-3/2}), \end{aligned} \quad (58)$$

with $\bar{Y} = 4b_0Y$ and $\xi_1 \approx -2.338$ is the rightmost zero of the Airy function. Setting $B_g = 0$ in this expression, we recover the DLA result obtained in [51]. Note also that, since $B_g = -1/b_0$, the coefficient of the $1/\bar{Y}$ term is equal to $-b_0(3 + 8b_0)$. Surprisingly, the $\mathcal{O}(Y^{-7/6})$ and $\mathcal{O}(Y^{-4/3})$ terms are not affected by the single log corrections (they do not depend on $B_g$). The NLL DGLAP contribution (which appears as a $1/\gamma$ pole in the NLL BFKL kernel [81]) starts contributing at order $Y^{-3/2}$, but this order is nonuniversal. The $2\zeta(3)\gamma^2$ term in $\chi_{LL}$ starts contributing at order $n = 6$ in the leading edge expansion, meaning that the correction of order $Y^{-2}$ in $\dot{\rho}_s$ would depend on $\zeta(3)$. The natural suppression of the NLL term in Eq. (51) by $1/Y^2$ due to the running of the coupling typically at the scale $Q_s$ implies that the universal terms of the asymptotic expansion of the saturation scale at large $Y$ only depends on $b_0$. In turn, higher loop corrections to the BFKL equation, such as $N^2$LL corrections which would start at order $\bar{\alpha}_s^3$, would be even more suppressed at large $Y$, so that they do not modify the universal series Eq. (58). A more formal argument illustrating the interplay between the $\alpha_s$ and $Y$ expansions goes as follows: including the $N^p$LL BFKL kernel at order $\alpha_s^{p+1}$ in the right-hand side of Eq. (51) gives the minimal $1/Y$ power arising from this contribution, which reads

$$\alpha_s^{p+1}(\rho_s(Y)) \tilde{\chi}_{N^p LL}(\beta/Y^{1/2}) \sim Y^{-(p+1)/2}, \quad (59)$$

since $\alpha_s(\rho_s(Y)) \sim 1/Y$ and $\tilde{\chi}_{N^p LL}(\gamma) \sim 1/\gamma^{p+1}$ at small $\gamma$. Thus, for $p \ge 2$, the $N^p$LL BFKL kernel contributes at most to the nonuniversal coefficients $\mathcal{O}(Y^{-(p+1)/2})$.

It is also interesting to compare with the universal asymptotic expansion of the saturation scale in the case of small-$x$ evolution (see [65,66,72,94,95]). For the energy dependence of $Q_s$, the universal terms also depend only on the LL and NLL BFKL evolution, but the coefficients of the development have a much stronger dependence on the shape of the LL and NLL BFKL kernel since they depend on the first five derivatives of the kernel at the small-$x$ saddle point $\beta_c = 0.6275$. Again, this is a specificity of the jet quenching problem, which is controlled by double logarithmic physics.

To sum up, the effect of the running coupling is to reduce the sensitivity of the subasymptotic corrections to higher orders in the resummation. Eventually, the universal asymptotic expansion of the saturation momentum is entirely given by the leading poles of the LL + NLL BFKL kernel in Mellin space or the singular plus finite part of the LO DGLAP splitting function. This is a consequence of the double logarithmic nature of the problem at hand.

A comparison between the double and single log result is shown Fig. 4 for two truncations of the asymptotic expansion, either up to the $1/Y$ term or up to the $\ln(Y)/Y^{3/2}$. In the latter case, one observes that the effect of single log corrections is very mild. The band show the estimated uncertainty coming from the following subleading correction in the large $Y$ development of $\dot{\rho}_s$. As expected, pushing the series up to order $1/Y^{3/2}$ reduces this source of uncertainty, and the corrections from single-log effects fall within the bands. However, we also know that the asymptotic series converges very slowly at small $Y \lesssim 10$ when using the nonlinear saturation boundary $\rho_s(Y)$ [51], so our estimation of the uncertainty coming

<figure><img src="image_6.png" /><figcaption>FIG. 4. The front wave velocity as a function of $Y$ in QCD. The bands show the uncertainty related to the $\bar{Y}^{-7/6}$ (dotted) and $\bar{Y}^{-3/2}$ (full) corrections in $\dot{\rho}_s$, obtained by adding the corresponding power multiplied by a coefficient $\kappa b_0$ with $|\kappa| \le 10$. This comprehensive interval for $\kappa$ is due to the potentially large nonuniversal $Y^{-3/2}$ coefficient, as noted in [51].</figcaption></figure>
from the nonuniversal $Y^{-3/2}$ term should be taken with a grain a salt in that domain.

## VI. SUMMARY AND OUTLOOK

In this paper, we have computed for the first time the universal behavior—independent of the nonperturbative, tree-level physics—of the saturation momentum associated with transverse momentum broadening of high energy partons in QCD media beyond the double logarithmic approximation. Our study relies on two pillars (i) the universal terms of the asymptotic series of $Q_s(L)$ at large $L$ are essentially controlled by the linearized evolution equation for the dipole cross section. This is a consequence of the mathematical mapping between this evolution equation and equations describing the propagation of traveling wave fronts into unstable states (ii) the quantum evolution of the dipole cross section in the jet quenching problem is dominated by the double logarithmic regime of pQCD since for $k_{\perp}^2 \sim Q_s^2$, $\ln(k_{\perp}^2/\mu^2) \sim \ln(L/\tau_0) \sim \ln(1/x)$. As a consequence, the universal asymptotic series can be obtained either from a BFKL or DGLAP approach, and the coefficient of the series are only sensitive to the shape of these two kernels close to the double logarithmic regime.

For a fixed coupling evolution, as in the case of the supersymmetric $N=4$ SYM theory, we have obtained the coefficients of the development of $\dot{\rho}_s$ up to 3-loops order $\alpha_s^{3/2}$, thanks to the known pole structure of BFKL at 3-loops in this theory. Using DGLAP/BFKL duality, we believe that our method can be straightforwardly extended to higher orders. As we observe a good convergence of the coefficients for $\alpha_s$ up to ~0.4, we have not tried to extend further our calculation.

The running coupling dramatically changes both the behavior of the large $L$ expansion and the way higher order corrections in the resummation appear. Interestingly enough, the universal terms depend only on the 1-loop QCD $\beta$ function. The main result of this paper, the asymptotic behavior of $d\rho_s/dY$ in QCD given by Eq. (58), is therefore exact to all orders in perturbation theory. At large $L$, the effect of corrections beyond the double logarithmic approximation turns out to be very mild, demonstrating the excellent convergence of $\rho_s$ after quantum evolution.

The main limitation of our work is the large $L$ assumption. Despite the calculation of the subasymptotic corrections which, in principle, enable us to reach phenomenological values for the system size $L$, our work should be supplemented by a study of the moderate and small $Y$ domain, which is not anymore driven by universality. This requires us to solve numerically a nonlinear evolution for the jet quenching parameter $\hat{q}$, valid at single logarithmic accuracy, along the lines of [42]. This is a very challenging task in practice, that we leave for future works.

Another possibility would be to take advantage of recent progress in quantum computing to address this problem, following the approach of [96,97], in which radiative corrections are straightforward to include.

## ACKNOWLEDGMENTS

We thank Edmond Iancu for insightful discussions about the strong coupling limit of the saturation scale. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-SC0012704. Y. M.-T. acknowledges support from the RHIC Physics Fellow Program of the RIKEN BNL Research Center.

## APPENDIX A: SINGLE LOGARITHMIC CORRECTIONS FROM DGLAP-LIKE EVOLUTION

In this appendix, we compute the single log corrections from a DGLAP evolution. Our starting point is the DGLAP equation for the gluon distribution function $xg(x, Q^2)$:

$$ \frac{\partial xg(x, Q^2)}{\partial \ln(Q^2)} = \frac{\alpha_s(Q^2)}{2\pi} \int_x^1 \frac{dz}{z} P_g(z) xg(x/z, Q^2). \quad (A1) $$

In the double logarithmic approximation, one uses $P_g(z) = 2C_A/z$, giving

$$ \frac{\partial xg(x, Q^2)}{\partial \ln(Q^2)} = \bar{\alpha}_s \int_x^1 \frac{dz}{z} \frac{x}{z} g(x/z, Q^2). \quad (A2) $$

To obtain the corresponding equation for $\hat{q}$, we use the fact that in the dilute limit, $\hat{q}(\rho, Y) \propto xg(x, Q^2)$ [51], with the identification $Y = -\ln(x)$ and $\rho = \ln(Q^2/\mu^2)$, so that the DGLAP evolution of $\hat{q}$ reads

$$ \frac{\partial \hat{q}(Y, \rho)}{\partial \rho} = \bar{\alpha}_s(\rho) \int_0^Y dY' \hat{q}(Y', \rho), \quad (A3) $$

which is almost equivalent to the nonlinear equation studied in [45] modulo the replacement $Y \to \min(Y, Y_s(\rho))$ in the upper limit of the $Y'$ integral [with $Y_s(\rho)$ the inverse function of $\rho_s(Y)$]. At single log accuracy, one can use the following approximation of $P_g(z)$:

$$ P_g(z) = \frac{2C_A}{z} (1 + B_g z), \quad (A4) $$

with $B_g = -11/12 - N_f/(6N_c^3)$ the finite part of the $g \to gg$ (and $g \to q\bar{q}$) splitting function. Hence, the equation we shall study is

$$ \frac{\partial \hat{q}(Y, \rho)}{\partial \rho} = \bar{\alpha}_s(\rho) \int_0^Y dY' [1 + B_g e^{Y'-Y}] \hat{q}(Y', \rho). \quad (A5) $$
The equation above can be written in a fully differential form:

$$ \frac{\partial^3 \hat{q}(Y, \rho)}{\partial Y^2 \partial \rho} + \frac{\partial^2 \hat{q}(Y, \rho)}{\partial Y \partial \rho} - \bar{\alpha}_s(\rho)(1 + B_g) \frac{\partial \hat{q}(Y, \rho)}{\partial Y} - \bar{\alpha}_s(\rho) \hat{q}(Y, \rho) = 0. \quad (A6) $$

Let us consider the fixed coupling approximation $\bar{\alpha}_s(\rho) = \bar{\alpha}_s$, and look for an asymptotic solution of the form

$$ \hat{q}(Y, \rho) = e^{\rho_s(Y)-Y} e^{\beta x}, \quad x = \rho - \rho_s(Y), \quad (A7) $$

with $\dot{\rho}_s = c$. We have

$$ \frac{\partial^2 \hat{q}(Y, \rho)}{\partial Y \partial \rho} = e^{\beta x} [(c-1)\beta - c\beta^2]. \quad (A8) $$

Plugging our ansatz inside the differential equation, one finds the following relation between $\beta$ and $c$:

$$ (\beta - 1)(1 + c(\beta - 1))\beta c + \bar{\alpha}_s(B_g + (\beta - 1)(1 + B_g)c) = 0. \quad (A9) $$

The extremum condition, obtained by differentiating the relation above with respect to $\beta$ gives

$$ -1 + \bar{\alpha}_s(1 + B_g) + 2\beta + c - 4\beta c + 3\beta^2 c = 0. \quad (A10) $$

This system can be exactly solved, but it is more enlightening to find the $\bar{\alpha}_s$ expansion of $c$ and $\beta$. The following two developments are solution

$$ c = 1 + 2\sqrt{\bar{\alpha}_s} + (2 + B_g)\bar{\alpha}_s + \mathcal{O}(\bar{\alpha}_s^{3/2}), \quad (A11) $$

$$ \beta = \sqrt{\bar{\alpha}_s} + (B_g - 1)\bar{\alpha}_s + \mathcal{O}(\bar{\alpha}_s^{3/2}). \quad (A12) $$

This is the same result as the one obtained from the fixed coupling NLL BFKL evolution.

## APPENDIX B: LEADING EDGE EXPANSION CALCULATION

We briefly reproduce here the calculation of the universal terms in the leading edge expansion detailed in [51]. We essentially detail the calculation of the first term $G_{-1}(z)$ in this series, which fixes the first subasymptotic correction to $\dot{\rho}_s$. The computation of the higher order terms proceeds in a similar fashion.

*Fixed coupling.*—We first insert Eq. (45) in Eq. (29), with the asymptotic series

$$ \dot{\rho}_s = c + \frac{\delta_1}{Y} + \dots \quad (B1) $$

Gathering the leading terms proportional to $1/Y$ in the resulting equation, we get the following differential equation satisfied by the function $G_{-1}$:

$$ -\frac{1}{2}\chi''(\beta_c)G_{-1}'' - \frac{1}{2}zG_{-1}' + \left(\frac{1}{2} + \delta_1 - \delta_1\beta_c\right)G_{-1} = 0. \quad (B2) $$

This equation can be solved analytically with the initial conditions $G_{-1}(z) = \text{uste} \times z + \mathcal{O}(z^2)$. This initial condition comes from the saturation constraint on the evolution equation which imposes $f(0, Y) = 1$ and therefore $G_{-1}(0) = 0$. Demanding the solution to decay at large $z$ and to be positive, one can further constrain the value of $\delta_1$ to be

$$ \delta_1 = \frac{3}{2(1 - \beta_c)}, \quad (B3) $$

so that the function $G_{-1}$ reads

$$ G_{-1}(z) = \text{uste} \times \beta_c z \exp \left( - \frac{z^2}{2\chi''(\beta_c)} \right). \quad (B4) $$

*Running coupling.*—The running coupling case is very similar, the only difference comes from the evolution equation (51) which imposes a diffusive power $1/6$ in the leading edge expansion. The homogeneity of Eq. (51) constrains the possible power of the subleading asymptotic corrections to $\dot{\rho}_s$. We have in particular,

$$ \dot{\rho}_s = 1 + \frac{c}{Y^{1/2}} + \frac{\delta_1}{Y^{5/6}} + \dots \quad (B5) $$

Plugging the series (57) and this development for $\dot{\rho}_s$ inside Eq. (51) and expanding in powers of $Y$ the result, we find the equation

$$ -G_{-1}''(z) + \left(\frac{1}{2}b_0z + \delta_1\sqrt{b_0}\right)G_{-1}(z) = 0, \quad (B6) $$

by identifying the $Y^{-2/3}$ power. The solution to this equation with boundary conditions $G_{-1}(z) \propto z$ at small $z$ and $G_{-1}(z) \to 0$ at large $z$ is (we denote $\text{Ai}$ the Airy function of the first kind)

$$ G_{-1}(z) = \text{uste} \times \text{Ai}(\xi_1 + 2^{-1/3}b_0^{1/3}z), \quad (B7) $$

with the constant $\delta_1$ fixed by these boundary conditions to

$$ \delta_1 = 2^{-2/3} b_0^{1/6} \xi_1. \quad (B8) $$

## APPENDIX C: COEFFICIENTS OF $\dot{\rho}_s$ IN PLANAR $\mathcal{N}=4$ SYM THEORY

In this appendix, we compute the $\alpha_s$ expansion of the coefficients in the development of $\dot{\rho}_s$ given by Eq. (46):

$$ \dot{\rho}_s(Y) = c + \frac{\delta_1}{Y} + \frac{\delta_2}{Y^{3/2}} + \mathcal{O}\left(\frac{1}{Y^2}\right), \quad (C1) $$
in the planar limit of the conformal $\mathcal{N} = 4$ SYM theory. At order $\alpha_s^{3/2}$, we have from the identities Eqs. (38) and (39),

$$c = 1 + 2\sqrt{\bar{\alpha}_s} + 2\bar{\alpha}_s + \left(1 - \frac{\pi^2}{6}\right)\bar{\alpha}_s^{3/2} + \mathcal{O}(\bar{\alpha}_s^2), \quad (C2)$$

$$\beta_c = \sqrt{\bar{\alpha}_s} - \bar{\alpha}_s + \left(\frac{1}{2} - \frac{\pi^2}{4}\right) \bar{\alpha}_s^{3/2} + \mathcal{O}(\bar{\alpha}_s^2). \quad (C3)$$

Therefore, the coefficient of the $1/Y$ term in $\dot{\rho}_s$ reads

$$\delta_1 = \frac{3}{2(\beta_c - 1)}, \quad (C4)$$

$$= -\frac{3}{2} - \frac{3}{2}\sqrt{\bar{\alpha}_s} + \left(\frac{3}{4} + \frac{3\pi^2}{8}\right)\bar{\alpha}_s^{3/2} + \mathcal{O}(\bar{\alpha}_s^2), \quad (C5)$$

and the coefficient of the $1/Y^{3/2}$ term is given by

$$\delta_2 = \frac{3}{2(\beta_c - 1)^2} \sqrt{\frac{2\pi}{\chi''(\beta_c)}}, \quad (C6)$$

$$= \frac{3\sqrt{\pi}}{2\bar{\alpha}_s^{1/4}} \left[ \bar{\alpha}_s^{1/2} + \frac{1}{2}\bar{\alpha}_s - \left(\frac{21}{8} - \frac{\pi^2}{8}\right)\bar{\alpha}_s^{3/2} + \mathcal{O}(\bar{\alpha}_s^2) \right]. \quad (C7)$$

[1] K. Adcox et al. (PHENIX Collaboration), Phys. Rev. Lett. 88, 022301 (2002).

[2] S. S. Adler et al. (PHENIX Collaboration), Phys. Rev. Lett. 91, 072301 (2003).

[3] C. Adler et al. (STAR Collaboration), Phys. Rev. Lett. 89, 202301 (2002).

[4] C. Adler et al. (STAR Collaboration), Phys. Rev. Lett. 90, 082302 (2003).

[5] J. Adams et al. (STAR Collaboration), Phys. Rev. Lett. 91, 172302 (2003).

[6] B. Abelev et al. (ALICE Collaboration), Phys. Lett. B 720, 52 (2013).

[7] B. Abelev et al. (ALICE Collaboration), J. High Energy Phys. 03 (2014) 013.

[8] S. Chatrchyan et al. (CMS Collaboration), Eur. Phys. J. C 72, 1945 (2012).

[9] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 114, 072302 (2015).

[10] J.-P. Blaizot and Y. Mehtar-Tani, Int. J. Mod. Phys. E 24, 1530012 (2015).

[11] G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E 24, 1530014 (2015).

[12] L. Cunqueiro and A. M. Sickles, Prog. Part. Nucl. Phys. 124, 103940 (2022).

[13] A. Adare et al. (PHENIX Collaboration), arXiv:1501.06197.

[14] Z. Citron et al., CERN Yellow Rep. Monogr. 7, 1159 (2019).

[15] J.-P. Blaizot, E. Iancu, and Y. Mehtar-Tani, Phys. Rev. Lett. 111, 052001 (2013).

[16] M. Aaboud et al. (ATLAS Collaboration), Phys. Lett. B 790, 108 (2019).

[17] S. Acharya et al. (ALICE Collaboration), Phys. Rev. C 101, 034911 (2020).

[18] P. Caucal, E. Iancu, and G. Soyez, J. High Energy Phys. 10 (2019) 273.

[19] Y. Mehtar-Tani, D. Pablos, and K. Tywoniuk, Phys. Rev. Lett. 127, 252301 (2021).

[20] A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 77, 011901 (2008).

[21] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 105, 252303 (2010).

[22] J. Adam et al. (ALICE Collaboration), J. High Energy Phys. 09 (2015) 170.

[23] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. C 96, 024905 (2017).

[24] A. H. Mueller, B. Wu, B.-W. Xiao, and F. Yuan, Phys. Lett. B 763, 208 (2016).

[25] L. Chen, G.-Y. Qin, S.-Y. Wei, B.-W. Xiao, and H.-Z. Zhang, Phys. Lett. B 773, 672 (2017).

[26] F. Ringer, B.-W. Xiao, and F. Yuan, Phys. Lett. B 808, 135634 (2020).

[27] M. S. Abdallah et al. (STAR Collaboration), Phys. Rev. C 105, 044906 (2022).

[28] F. D'Eramo, M. Lekaveckas, H. Liu, and K. Rajagopal, J. High Energy Phys. 05 (2013) 031.

[29] F. D'Eramo, K. Rajagopal, and Y. Yin, J. High Energy Phys. 01 (2019) 172.

[30] G. Moliere, Z. Naturforsch. A 3, 78 (1948).

[31] J. a. Barata, Y. Mehtar-Tani, A. Soto-Ontoso, and K. Tywoniuk, Phys. Rev. D 104, 054047 (2021).

[32] A. V. Sadofyev, M. D. Sievert, and I. Vitev, Phys. Rev. D 104, 094044 (2021).

[33] J. a. Barata, A. V. Sadofyev, and C. A. Salgado, Phys. Rev. D 105, 114010 (2022).

[34] C. Andres, F. Dominguez, A. V. Sadofyev, and C. A. Salgado, Phys. Rev. D 106, 074023 (2022).

[35] G. D. Moore and N. Schlusser, Phys. Rev. D 101, 014505 (2020); 101, 059903(E) (2020).

[36] G. D. Moore, S. Schlichting, N. Schlusser, and I. Soudi, J. High Energy Phys. 10 (2021) 059.

[37] T. Liou, A. H. Mueller, and B. Wu, Nucl. Phys. A916, 102 (2013).

[38] J.-P. Blaizot, F. Dominguez, E. Iancu, and Y. Mehtar-Tani, J. High Energy Phys. 06 (2014) 075.

[39] S. Caron-Huot, Phys. Rev. D 79, 065039 (2009).
[40] J. Ghiglieri and E. Weitz, J. High Energy Phys. 11 (2022) 068.

[41] J.-P. Blaizot and Y. Mehtar-Tani, Nucl. Phys. A929, 202 (2014).

[42] E. Iancu, J. High Energy Phys. 10 (2014) 095.

[43] P. Arnold, J. High Energy Phys. 03 (2022) 134.

[44] P. Arnold, T. Gorda, and S. Iqbal, J. High Energy Phys. 04 (2022) 085.

[45] P. Caucal and Y. Mehtar-Tani, Phys. Rev. D 106, L051501 (2022).

[46] G. Dee and J. Langer, Phys. Rev. Lett. 50, 383 (1983).

[47] W. Van Saarloos, Phys. Rev. Lett. 58, 2571 (1987).

[48] S. Munier and R. B. Peschanski, Phys. Rev. Lett. 91, 232001 (2003).

[49] U. Ebert and W. van Saarloos, Physica (Amsterdam) 146D, 1 (2000).

[50] W. Van Saarloos, Phys. Rep. 386, 29 (2003).

[51] P. Caucal and Y. Mehtar-Tani, J. High Energy Phys. 09 (2022) 023.

[52] E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, Sov. Phys. JETP 45, 199 (1977).

[53] I. I. Balitsky and L. N. Lipatov, Sov. J. Nucl. Phys. 28, 822 (1978).

[54] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 781 (1972) [Sov. J. Nucl. Phys. 15, 438 (1972)].

[55] G. Altarelli and G. Parisi, Nucl. Phys. B126, 298 (1977).

[56] Y. L. Dokshitzer, Zh. Eksp. Teor. Fiz. 73, 1216 (1977) [Sov. Phys. JETP 46, 641 (1977)].

[57] E. Iancu and D. N. Triantafyllopoulos, Phys. Rev. D 90, 074002 (2014).

[58] H. Kowalski and D. Teaney, Phys. Rev. D 68, 114005 (2003).

[59] T. Lappi, Phys. Lett. B 703, 325 (2011).

[60] P. B. Arnold and W. Xiao, Phys. Rev. D 78, 125008 (2008).

[61] A. Peshier, Phys. Rev. Lett. 97, 212301 (2006).

[62] W. A. Horowitz and Y. V. Kovchegov, Nucl. Phys. A849, 72 (2011).

[63] Y. V. Kovchegov and H. Weigert, Nucl. Phys. A807, 158 (2008).

[64] S. Munier, Phys. Rep. 473, 1 (2009).

[65] A. H. Mueller and D. N. Triantafyllopoulos, Nucl. Phys. B640, 331 (2002).

[66] G. Beuf, arXiv:1008.0498.

[67] F. Dominguez, A. H. Mueller, S. Munier, and B.-W. Xiao, Phys. Lett. B 705, 106 (2011).

[68] F. Dominguez, J.-W. Qiu, B.-W. Xiao, and F. Yuan, Phys. Rev. D 85, 045003 (2012).

[69] V. Vaidya, arXiv:2107.00029.

[70] V. Vaidya, arXiv:2109.11568.

[71] V. S. Fadin and L. N. Lipatov, Phys. Lett. B 429, 127 (1998).

[72] S. Munier and R. B. Peschanski, Phys. Rev. D 69, 034008 (2004).

[73] Y. V. Kovchegov and E. Levin, *Quantum Chromodynamics at High Energy*, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University Press, Cambridge, England, 2012).

[74] G. P. Salam, J. High Energy Phys. 07 (1998) 019.

[75] M. Ciafaloni, D. Colferai, and G. P. Salam, Phys. Rev. D 60, 114036 (1999).

[76] G. Altarelli, R. D. Ball, and S. Forte, Nucl. Phys. B742, 1 (2006).

[77] G. Beuf, Phys. Rev. D 89, 074039 (2014).

[78] B. Ducloué, E. Iancu, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, J. High Energy Phys. 04 (2019) 081.

[79] I. Balitsky, Nucl. Phys. B463, 99 (1996).

[80] Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999).

[81] S. Marzani, R. D. Ball, P. Falgari, and S. Forte, Nucl. Phys. B783, 143 (2007).

[82] M. S. Costa, V. Goncalves, and J. Penedones, J. High Energy Phys. 12 (2012) 091.

[83] N. Gromov, F. Levkovich-Maslyuk, and G. Sizov, Phys. Rev. Lett. 115, 251601 (2015).

[84] V. N. Velizhanin, arXiv:1508.02857.

[85] S. Caron-Huot and M. Herranen, J. High Energy Phys. 02 (2018) 058.

[86] J. Ghiglieri and H. Kim, J. High Energy Phys. 12 (2018) 049.

[87] A. V. Kotikov and L. N. Lipatov, Nucl. Phys. B582, 19 (2000).

[88] Y. Hatta, E. Iancu, and A. H. Mueller, J. High Energy Phys. 01 (2008) 063.

[89] Y. Hatta, E. Iancu, and A. H. Mueller, J. High Energy Phys. 01 (2008) 026.

[90] F. Dominguez, C. Marquet, A. H. Mueller, B. Wu, and B.-W. Xiao, Nucl. Phys. A811, 197 (2008).

[91] V. N. Velizhanin, arXiv:2106.06527.

[92] A. H. Mueller, B. Wu, B.-W. Xiao, and F. Yuan, Phys. Rev. D 95, 034007 (2017).

[93] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997).

[94] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rep. 100, 1 (1983).

[95] E. Iancu, K. Itakura, and L. McLerran, Nucl. Phys. A708, 327 (2002).

[96] M. Li, T. Lappi, and X. Zhao, Phys. Rev. D 104, 056014 (2021).

[97] J. a. Barata, X. Du, M. Li, W. Qian, and C. A. Salgado, Phys. Rev. D 106, 074013 (2022).