BFKL and CCFM evolutions with saturation boundary

We perform numerical studies of the BFKL and CCFM equations for the unintegrated gluon distribution supplemented with an absorptive boundary which mimics saturation. For the BFKL equation, this procedure yields the same results for the saturation momentum and the gluon distribution above saturation as the non-linear BK equation, for both fixed and running coupling, and for all the considered energies. This similarity goes beyond expectations based on the correspondence with statistical physics, which hold only for fixed coupling and asymptotically high energies. For the CCFM equation, whose non-linear generalization is not known, our method provides the first study of the approach towards saturation. We find that, in the running-coupling case, the CCFM and BFKL predictions for the energy dependence of the saturation momentum are identical within our numerical accuracy. A similar saturation boundary could be easily implemented in the CCFM-based Monte Carlo event generators, so like CASCADE.


Introduction
The imminent high-energy experiments at LHC will considerably enlarge the phase-space where the unitarity corrections to QCD interactions, like gluon saturation and multiple scattering, are expected to be important. Such corrections should in particular influence some 'hard' observables, like jet production at forward rapidities, whose theoretical description lies within the realm of perturbative QCD. The jets to be measured at LHC will carry relatively large transverse momenta Q ≥ 10 GeV, but because of the high-energy kinematics, their description may go beyond the standard pQCD formalism at high Q 2 -the DGLAP evolution [1] of the parton distributions together with the collinear factorization of the hadronic cross-sections. Rather, the high-energy evolution, of the BFKL [2] or CCFM [3] type, and the associated k T -factorization scheme should prevail whenever the energy logarithms ln s ∼ Y are larger than the momentum ones, ln Q 2 . Besides, this evolution is expected to be amended by non-linear effects reflecting gluon saturation and multiple scattering, whose theoretical description within pQCD has been given only recently [4,5,6,7,8] and (in full rigor) only in the leading-order approximation (see Refs. [9,10] for recent extensions to running coupling). Such effects can make themselves felt even at relatively large momenta Q, well above the saturation momentum Q s (the characteristic scale for the onset of unitarity corrections), via phenomena like geometric scaling [11,12,13,14,15,16,17], which reflect the change in the unintegrated gluon distribution at high k ⊥ ≫ Q s due to saturation at low k ⊥ Q s . The saturation scale Q s grows, roughly, like a power of the energy: Q 2 s ∼ s λ with λ ≃ 0.25 from fits to HERA data [18,19], which are also supported by next-to-leading order theoretical calculations [14]. For forward jet production in proton-proton collisions at LHC, Q s is expected in the ballpark of 2 to 3 GeV. Besides, much higher values of Q s can be effectively reached [20] by focusing on 'hot spots' (partons at large x and high Q 2 , which develop their own gluon cascades and saturation momentum) in some rare events, so like Mueller-Navelet jets.
In view of the above, it becomes important and urgent to provide realistic, quantitative, predictions for the effects of saturation on relatively hard (Q 2 ≫ Q 2 s ) observables at LHC. Besides making the whole perturbative approach fully justified, the restriction to relatively hard momenta entails some important simplifications, which are essential for the strategy that we shall propose in this Letter. First, this implies that one can neglect the complex many-body correlations which develop at saturation (Q Q s ), but focus on the gluon phase-space density (or 'unintegrated gluon distribution') alone. Hence, the standard k T -factorization of the hadronic cross-sections still applies, but with modified unintegrated gluon distributions, which reflect saturation. This opens the possibility to include the effects of saturation within Monte-Carlo event generators based on k T -factorization, so like CASCADE [21], which relies on the CCFM evolution.
Second, this means that the saturation effects in the gluon distribution and, in particular, the saturation momentum itself can be computed without a detailed knowledge of the non-linear dynamics responsible for unitarization. Rather, they are fully determined by the linear evolution if the latter is supplemented with an absorptive boundary condition at low momenta, whose position is energydependent (as this mimics the saturation momentum) and is self-consistently determined by the evolution [12,13,14] . This property is both important and highly non-trivial. It is important as it allows one to perform studies of saturation even for evolution equations whose non-linear generalizations are not known, so like the CCFM evolution, and also the BFKL evolution beyond the leading-order approximation. It is moreover non-trivial since the high-energy evolution is non-local in transverse momenta, hence the growth in the gluon distribution at high momenta k ⊥ ≫ Q s could be well feed by radiation from lower momenta k ⊥ Q s . This is clearly the case for the linear evolution with running coupling, in which the gluon distribution grows faster in the infrared (where the coupling is stronger) and then acts as a source for radiating gluons at higher momenta.
Yet, at least for a fixed coupling and for asymptotically high energies, it has been firmly established that the high-energy evolution towards saturation is driven by the linear evolution in the dilute tail of the gluon distribution at high transverse momenta k ⊥ ≫ Q s . The respective argument is based on a correspondence between high-energy QCD and statistical physics [15,22], which however fails to apply in the more realistic case of a running coupling [23]. For that case, our subsequent numerical results will provide the first unambiguous evidence that the evolution of the saturation front in QCD is indeed driven by the linear part of the evolution, at all energies, including the preasymptotic ones. Most precisely, we shall find that the BFKL equation with saturation boundary provides exactly the same results for the saturation momentum Q s (Y ) and for the gluon distribution at k ⊥ ≥ Q s (Y ) as the non-linear Balitsky-Kovchegov (BK) equation [4,5] , for both fixed and running coupling, and for all the considered rapidities Y ≤ 120 -including lower values Y ≤ 14, as relevant for the phenomenology at LHC.
But the BFKL equation and its non-linear, BK, extension, to be discussed in Sect. 2, will merely serve as a playground to test our method for numerically implementing the saturation boundary condition within a generic linear evolution. Our main interest is rather in the CCFM evolution, that we shall discuss in Sect. 3, and which stays at the basis of Monte-Carlo event generators [21]. There are at least two reasons why the CCFM evolution is a privileged tool in that respect. First, it takes into account the quantum coherence between successive emissions, leading to angular ordering in the parton cascades. This allows for a more realistic description of the particle distribution in the final state as compared to the BFKL evolution. (The latter only guarantees the correct treatment of inclusive quantities like the unintegrated gluon distribution.) Second, the CCFM equation provides an interpolation between BFKL dynamics at small x and (approximate) DGLAP dynamics at larger x. Of course, when discussing saturation we are a priori interested in the small-x region, and there is little doubt that, for asymptotically high energies, the approach towards saturation is more correctly described by the BFKL formalism. However, the energies to be available at LHC are far from asymptotia, and a better treatment of the transverse momentum ordering within the parton cascades, as implicit in the CCFM formalism, is probably essential for most of the "small-x" phase-space to be experimentally accessible at LHC.
In this Letter, we shall limit ourselves to applications of the CCFM and BFKL evolutions of the unintegrated gluon distribution, so it will be meaningful to compare their results. To achieve a faster numerical convergence, we shall solve the respective differential equations, rather than building Monte-Carlo codes. But our prescription for implementing the saturation boundary condition is straightforward to apply within a Monte-Carlo event generator, so like CASCADE. Our present analysis can be viewed as a first step towards building an event generator which includes saturation [24]. The CCFM equation is considerably more involved than the BFKL one (in particular, its solution involves one additional variable: a maximum angle), so in order to be able to follow this evolution up to relatively high rapidities we shall consider a slightly simplified version of it. The simplification is obtained by rewriting the CCFM equation in a 'more inclusive form', i.e., by using the virtual, Sudakov-like, terms to cancel some of the real gluon emissions [25,26]. This procedure also involves some kinematical approximations, which however are in the spirit of the CCFM formalism 1 . By solving the ensuing equation with saturation boundary condition, we shall for the first time study the onset of unitarity corrections within the CCFM evolution. One of our most interesting results is that, in the running coupling case, the CCFM evolution in the presence of saturation provides almost identical results (for the saturation momentum and the gluon spectrum above Q s ) as the respective BFKL evolution.

BFKL evolution with absorptive boundary
In this section we shall explain our method for effectively implementing saturation within a unitarityviolating linear evolution on the example of the BFKL equation [2]. This is interesting since the corresponding non-linear equation which obeys unitarity is known as well, the BK equation [4,5], and thus it can be used to test our method. Although all our numerical studies will be performed in (transverse) momentum space (referring to the unintegrated gluon distribution), it is more convenient to explain our method in coordinate space. Then, the BK equation describes the high-energy evolution of the scattering amplitude T (Y, r) of a small quark-antiquark dipole with transverse size r. We shall assume the target to be infinite and homogeneous in transverse directions, so we can ignore the impact-parameter dependence of the scattering amplitude and average over angles. The corresponding equation reads Hereᾱ s ≡ α s N c /π, and z and r − z are the transverse sizes of the two dipoles into which the parent dipole r has dissociated before scattering off the target. The last term, quadratic in T , in the r.h.s. of the equation describes multiple scattering (the simultaneous scattering of both daughter dipoles) and is responsible for unitarization. With this last term omitted, (2.1) reduces to the BFKL equation, which describes the unlimited (exponential) growth of the scattering amplitude with Y and the symmetric expansion of the support of T (Y, r) in r towards both small and large dipole sizes. Note however that the transverse non-locality in Eq. (2.1) is quite weak -the large daughter dipoles with size z ≫ r are suppressed by the 'dipole kernel' r 2 /[z 2 (r − z) 2 ], whereas the very small ones, with z ≪ r, are disfavoured by the T -dependent terms in the r.h.s., which exactly cancel each other as z → 0 -and can be described as diffusion in the logarithmic variable ρ ≡ ln(r 2 0 /r 2 ). Here, r 0 is an arbitrary scale of reference (say, the unitarization scale in the target at Y = 0).
However, the fully non-linear equation (2.1) preserves (and actually saturates) the unitarity bound T ≤ 1, as it manifestly has T = 1 as a fixed point. Because of that, the respective evolution is asymmetric in r (or ρ), and the solution T (Y, r) ≡ T (Y, ρ) looks like a front, which interpolates between T = 1 at relatively small ρ and T = 0 at ρ → ∞, and which with increasing Y propagates towards larger values of ρ. Behind the front, the scattering amplitude has reached the 'black disk' limit T = 1 and thus cannot grow anymore. Ahead of the front, the amplitude is still weak, T ≪ 1, so the non-linear term in Eq.
The previous discussion already suggests that the progression of the saturation front towards larger values of ρ is driven by the BFKL evolution of the dilute tail at ρ ≫ ρ s (Y ) -the front is 'pulled' by its tail. This is an important property, as it allows us to determine the position ρ s (Y ) of the front and its shape around ρ s (Y ) from studies of the linear, BFKL, evolution alone. This property is highly nontrivial, in view of the non-locality of the non-linear equation (2.1). Other non-linear equations which are non-local and exhibit saturation are known to develop a pushed front, i.e., a front whose progression is driven by the growth and accumulation of 'matter' behind the front [28].
For the case of the BK equation with fixed coupling and for sufficiently high energy, the pulledfront property follows from the identification, made in Ref. [15], between the asymptotic form of the BK equation at high energy 2 and the FKPP equation (from Fisher Kolmogorov, Petrovsky, and Piscounov) of statistical physics, for which this property has been established with mathematical rigor 3 . However, this identification does not extend to a running coupling, and hence it fails to apply for the real QCD problem.
The running of the coupling is known to have dramatic consequences for the high-energy evolution [12,13,14,23], and in particular for the approach towards saturation: the growth of the coupling with decreasing momenta, or increasing dipole sizes, amplifies the contribution of the latter to the evolution, which then becomes asymmetric even in the absence of saturation. In fact, the BFKL evolution with running coupling is infrared-unstable, in the sense that it requires an infrared cutoff to avoid the blowup of the QCD coupling at k ⊥ ∼ Λ QCD , and then the results of the evolution are strongly sensitive to the value of this cutoff -so that the whole procedure has no predictive power (see Fig. 3

below).
In that scenario (BFKL with running coupling), the growth of the gluon distribution at high momenta k ⊥ ≫ Λ QCD is mostly feeded by radiation from the bulk of the distribution at infrared (k ⊥ ∼ Λ QCD ) momenta. Hence, the whole perturbative framework becomes questionable and, besides, one may expect the associated saturation front -as generated after enforcing unitarity -to be of the 'pushed' type.
Yet, as our explicit numerical solutions will demonstrate, this is actually not the case: the saturation front remains of the 'pulled' type even for a running coupling. This is so because the gluon modes with k ⊥ Q s (Y ) become inert due to saturation, so the evolution is again driven by the dilute tail at high momenta, so like for fixed coupling. In particular, the infrared problem is automatically avoided: the saturation scale effectively acts as an infrared cutoff, which becomes 'hard' (Q 2 s (Y ) ≫ Λ 2 QCD ) for sufficiently high energy. This opens the way towards realistic studies of the front dynamics within the context of the linear evolution, as originally suggested in Refs. [12,13] . To that aim, the linear evolution equations must be supplemented with an appropriate saturation boundary condition, that we now describe.
Such a boundary condition must ensure that the amplitude never becomes bigger than one. By itself, the position of the front is not a priori known, but must be determined when solving the equation. To that aim, let us first introduce a line of constant amplitude ρ = ρ c (Y ) via the condition where the number c is strictly smaller than one, but not much smaller. (The saturation line ρ s (Y ) would correspond to c ∼ 1.) For ρ < ρ c (Y ) and sufficiently high energy, the solution T BFKL (Y, ρ) to the BFKL equation would become larger than one -in fact, arbitrarily large. If this equation is to be solved numerically, one may think about identifying the point ρ c (Y ) numerically at each step in Y , and then enforcing the unitarity limit T = 1 for any ρ which is smaller than ρ c (Y ) and sufficiently far away from it -say, for ρ ≤ ρ c (Y ) − ∆ with ∆ ≃ ln(1/c) a number of O(1). However, this would not be a very good strategy in practice, since T = 1 is not a fixed point for the BFKL equation, so an amplitude of O(1) would be exponentially amplified by the subsequent evolution. Even if, at small ρ, one cuts off this evolution by hand step-by-step, it is not clear (especially for running coupling) whether the spurious radiation from small ρ will not affect the tail of the front at large ρ. It is therefore preferable, as originally suggested in Ref. [13] , to enforce the amplitude to vanish for ρ ≤ ρ c (Y ) − ∆ : Since it yields, by construction, T = 0 beyond the saturation front, this procedure cannot be used for any physical problem which is sensitive to the black disk limit, like deep inelastic scattering at low Q 2 Q 2 s (Y ), or particle production at low transverse momenta. On the other hand, as we shall see, this procedure accurately describes the dynamics of the front, meaning its position and shape for ρ > ρ c (Y ), and hence it correctly provides the tail of the gluon distribution at transverse momenta , we see that our results cover the phenomenologically interesting region where geometric scaling is expected at high Q 2 (cf. Introduction).
To describe our numerical results, let us first change from the coordinate to the momentum representation, i.e. from the scattering amplitude T (Y, r) to the 'unintegrated gluon distribution' A(Y, k)the quantity which enters the calculation of cross-sections within the k T -factorization. For the present purposes, A(Y, k) can be defined as the following Fourier transform of the dipole amplitude [12] A(Y, k) = d 2 r 2πr 2 e −ik·r T (Y, r) .

(2.4)
With this definition, the more standard, 'integrated', gluon distribution is obtained as For our homogeneous target, A(Y, k, b) ≡ A(Y, k), hence the above integral over b simply yields the hadron transverse area πR 2 . With these conventions, the gluon occupation number -i.e., the number of gluons of a given color per unit rapidity per unit volume in transverse phase-space -is not exactly A(Y, k), but rather A(Y, k)/ᾱ s (up to a numerical factor). Via Eq. (2.4), the saturation front for T (Y, r) translates into a corresponding front for A(Y, k). For very large momenta, k ≫ Q s (Y ), one can use T ∼ r 2 ln(1/r 2 ) ('color transparency') and then Eq. (2.4) reproduces the bremsstrahlung spectrum: A(Y, k) ∼ 1/k 2 , or 4 A(Y, ρ) ∼ e −ρ . This behaviour is modified when k gets closer to (but still larger than) Q s , due to the BFKL evolution in the presence of saturation. For instance, in the fixed coupling case and for sufficiently high energy (ᾱ s Y ≫ 1), one finds [12,13,15] The exponential ∼ e −γsρ describes the modification of the bremsstrahlung spectrum due to the BFKL 'anomalous dimension' 1 − γ s , while the last, Gaussian, factor describes BFKL diffusion. The numbers γ s ≈ 0.63 and β ≈ 48.52 are specific for the problem at hand: they characterize the BFKL evolution along a line of constant amplitude, cf. Eq. (2.2) (see [12,13] for details). Also, the presence of the saturation scale ρ s and the overall factor ρ − ρ s in the r.h.s. of Eq. (2.6) are the hallmarks of saturation, and can 4 From now on, we shall use the notation ρ for either ln(r 2 0 /r 2 ), or ln(k 2 /k 2 0 ) (with k 0 = 1/r 0 ), the difference being clear from the context. be generated within the framework of the linear evolution only after imposing the saturation boundary condition (2.3) [13] . Under the same conditions, the saturation scale is obtained as where λ ≈ 4.88 and the last, constant, term is not under control (since Eq. (2.7) is merely an asymptotic expansion at high-energy).
In the opposite limit of very low momenta k ≪ Q s (Y ), the integral in Eq. (2.4) is dominated by large dipole sizes r ≫ 1/Q s for which T = 1; one then finds A(Y, k) ≃ ln[Q s (Y )/k] = [ρ s (Y ) − ρ]/2. Thus, behind the front, A(Y, k) is not exactly constant (unlike the dipole amplitude), but it is slowly growinglogarithmically in both 1/k and Y . What saturates at high density is not the gluon occupation number A/ᾱ s , but rather the rate for gluon emission [8].
The analytic results in Eqs. (2.6) and (2.7), together with the corresponding ones at running coupling [13,14] , have been already tested in the literature against numerical solutions to the BK equation. Here, we are rather interested to compare the solutions to the latter against numerical solutions for the BFKL equation supplemented with the boundary condition (2.3) -with both fixed and running coupling. The momentum-space version of the BK equation is particularly simple in that the non-linear term is local: The linear version of this equation, i.e., the BFKL equation for A(Y, k), will be solved with an absorptive boundary condition similar to that for T (Y, r) in Eq. (2.3). This is appropriate since, with the normalization in Eq. (2.4), the saturation effects in the gluon distribution become important when A(Y, k) ∼ O(1).
For the fixed coupling calculations we shall useᾱ s = 0.2. To include running coupling effects, we shall pull theᾱ s factor inside the q-integral in Eq. (2.8) and use the one-loop expression for the running coupling with scale Q 2 = max(k 2 , q 2 ) and Λ QCD = 200 MeV. This simple prescription is in agreement with the recently constructed running-coupling version of the BK equation [9,10] . To avoid the infrared divergence of the coupling at Q 2 = Λ 2 QCD , we shall replace α s (Q 2 ) → α s (Q 2 + µ 2 ) for some parameter µ. Our default choice will be µ 2 = 0.5 GeV 2 , but we shall study the sensitivity of our results to variations in µ. Our initial condition A(Y = 0, k) is given by the bremsstrahlung spectrum for k > 1 GeV (with maximal height A = 0.5) and it vanishes for k < 1 GeV.
For the fixed coupling case, our results are displayed in Fig. 1 for five values of the rapidity within the range 20 ≤ Y ≤ 70 and for two different choices for the parameters c and ∆. The most important observation about these results is that the saturation fronts generated by the two types of evolution do precisely coincide with each other for all momenta ρ ≥ ρ c (Y ) − ∆, and for all the considered rapidities. This property is not altered by changing the values for c and (correlated to it) for ∆. We have also checked that these numerical curves are consistent with the analytic estimates in Eqs. (2.6) and (2.7), including the expected values for λ and γ s . We now turn to the more realistic case of a running coupling. Then, as alluded to before, the pure BFKL evolution is infrared unstable, since the rapid growth of the gluon distribution at small values of k (where is coupling is larger) is also feeding the growth at higher k. Therefore the linear evolution behaves quite differently compared to the non-linear one, even at high k. This is clearly visible in Fig. 2 where we compare the strict BFKL evolution to the BK one, and to the BFKL evolution with the absorptive boundary. Once again, there is a perfect matching between the saturation fronts provided by BK and, respectively, BFKL with saturation boundary. (For the latter, we used the same values for c and ∆ as at fixed coupling.) On the other hand, we see a dramatic difference with respect to the linear evolution, which progresses much more rapidly towards the right.
From these curves, it is also possible to extract the Y -dependence of the saturation momentum ρ s (Y ) for running coupling. We find that the squared-root law ρ s ≃ λ r √ Y predicted by the theory [12,13,14] for asymptotically high energies provides a good fit to our numerical results for Y ≥ 10, with a fitted value λ r ≃ 2.8 which agrees reasonably well with the (asymptotic) theoretical expectation 5 λ r ≃ 3.2. Now, from the phenomenological point of view, we are more interested in values of Y which are not that large, say Y ≤ 14 (corresponding to x 10 −6 ), as relevant for forward jet production at LHC. With that in mind, we also show in Fig. 2 (right) the results for lower values of Y , between 6 and 14 units; one can thus see that the absorptive boundary method works equally well also for such lower rapidities.
Finally, to illustrate the infrared stability introduced by saturation, we exhibit in Fig. 3 results obtained for different values of the IR cutoff µ 2 inserted in the running coupling. Unlike the pure BFKL results (left figure), which are extremely sensitive to a change in µ, the results corresponding to the saturation boundary condition (right figure) show no sensitivity whatsoever. 5 Specifically, for asymptotically large Y , the running-coupling BFKL evolution yields [12,13] : ρs(Y ) ≃ √ 2λb 0 Y where λ ≃ 4.88 is the same number as in Eq. (2.7) and b 0 ≡ 12Nc/(11Nc − 2N f ) is the coefficient in the one-loop running coupling:ᾱ(Q 2 ) = b 0 /ln(Q 2 /Λ 2 QCD ). In our simulations, we use N f = 0, hence we expect λr ≡ √ 2λb 0 ≃ 3.26, which is indeed consistent with the fit to the curves in Fig. 2.

CCFM evolution with absorptive boundary
In this section we shall present a compact version of the CCFM equation [3] to which we shall apply the boundary condition described in the previous section. A more comprehensive discussion of the CCFM formalism and its relation to BFKL will be given elsewhere [27], together with more detailed numerical studies, of which the present Letter is only giving a glimpse.
As mentioned in the introduction, the CCFM evolution takes into account the quantum coherence between successive emissions via angular ordering in the parton cascades. Accordingly, the respective gluon distribution now depends on three variables, A = A(x, k,q), where the third variableq is a transverse momentum related to the maximum angle which determines the phase space where emissions are allowed. This angle is set by the hard scattering of the space-like photon against a quark inside the proton. It is customary to define the variable ξ which is the squared angle, ξ ≡ q 2 /(y 2 E 2 ), where q is the transverse momentum of a gluon emitted in the s-channel, y is its longitudinal momentum fraction, and E is the energy of the proton; then, all emissions must satisfy ξ ≤ξ ≡q 2 /(x 2 E 2 ).
The CCFM equation for A can be written in different versions, depending on how 'exclusive' we choose the gluon distribution to be 6 . That is, so long as one is not interested in the structure of the final state, one can 'integrate out' some of the emissions in the s-channel, as they do not change the overall (unintegrated) gluon distribution (and hence neither the probability for the interaction with a projectile). In practice, this amounts to suitable cancelations between 'real' gluon emissions and 'virtual' terms (or 'Sudakov factors'). Of course, if one is interested in studying the exclusive final states, all emissions which were removed from the initial state must later be included as final state radiation.
If one keeps within A only the emissions associated with the 1/z pole in the splitting function, then the CCFM equation can be written as the following integral equation where q and 1 − z are, respectively, the transverse momentum and the energy fraction of the 'real' gluon emitted (in one step of the evolution) in the s-channel 7 . The theta function comes from the angular ordering constraint ξ ≤ξ. There is also an energy ordering implicit in (3.9): emissions are ordered in where κ ≡ min(k 2 , k ′2 )/max(k 2 , k ′2 ) and Notice that h(κ) → 0 as κ → 1, so Eq. (3.14) has indeed no singularity at k ′ = k. To obtain an integrodifferential equation, we define z = e −y and x = e −Y and differentiate the l.h.s. with respect to Y . We thus get which is our final version for the CCFM equation and is equivalent to an equation originally proposed in Ref. [25] (although there this was not derived from the CCFM equation (3.9) in the way we did). The fact that this equation is nonlocal in Y should not come as a surprise, since the CCFM evolution is not really an evolution in Y , but rather in ξ.
It is first interesting to compare the predictions of Eq. (3.16) to those of the BFKL equation for the strictly linear evolution. This is shown in Fig. 4 for both fixed and running coupling, and one clearly sees  that the BFKL evolution is considerably faster. This difference is to be attributed to the non-local term in the r.h.s. of Eq. (3.16): the 'retarded' distribution A(Y − log(k ′2 /k 2 ), k ′ ) is generally smaller than the 'instantaneous' one A(Y, k ′ ).
But even though the CCFM evolution is somewhat slower, Eq. (3.16) still shows a pronounced growth with Y , which in the absence of any non-linearity would rapidly lead to unitarity violation. To cure for that, we now enforce the absorptive boundary condition (2.3) on Eq. (3.16). The corresponding results are compared to those of the purely (CCFM) linear evolution in Fig. 5, for both fixed and running coupling. As in the BFKL case (compare to Fig. 2), the difference is more pronounced for a running coupling, since then the linear evolution is again infrared-unstable, and this instability is removed by the inclusion of saturation.
For the more realistic, running-coupling, case it is furthermore interesting to show the results at lower rapidities, as relevant for LHC. This is exhibited in the leftmost figure in Fig. 6, together with the corresponding results of the strictly linear evolution. As one can see there, for Y = 14 and k as high as 10 GeV (which is well above the respective saturation momentum Q s ≃ 2.5 GeV), saturation reduces the predicted gluon distribution by about a factor of 2.
Finally, in the rightmost figure in Fig. 6, we compare the saturation fronts provided by the BFKL and CCFM evolutions with running coupling and for relatively high rapidities, up to Y = 120. What is remarkable about this last figure is that the BFKL and CCFM evolutions with saturation and running coupling appear to very close to each other and for all rapidities -meaning that the corresponding fronts propagate at roughly the same speed. This is at variance with the corresponding situation at fixed coupling, where the BFKL evolution is still faster. This is confirmed by the estimates of the saturation momentum for the CCFM evolution, as extracted from fits to our numerical results: It is numerically clear that ρ s rises linearly with Y for a fixed coupling, and it exhibits a √ Y behaviour for a running coupling, just like for BFKL. Trying simple fits of the form ρ s = λ fᾱs Y and, respectively, ρ s = λ r √ Y , we find the values 10 λ f ≈ 3.5 and, respectively, λ r ≈ 2.8. For fixed coupling, this value λ f for the saturation exponent is indeed smaller (although not much smaller) than the corresponding BFKL estimate λ f ≃ 4.9. But for running coupling, the CCFM and BFKL estimates for λ r are essentially the same within our numerical accuracy. Moreover, the manifest similarity between the shapes of the BFKL and CCFM fronts in Fig. 6 (right) suggests that the CCFM fronts exhibit the same properties of geometric scaling as the BFKL ones (again, for a running coupling).
It would be of course interesting to understand this similarity between the BFKL and CCFM evolutions towards saturation in more depth, and also to perform more detailed studies of the CCFM evolution with saturation boundary, in order e.g. to explicitly test geometric scaling. We postpone such studies to a further work [27].