Non-linear evolution in CCFM: The interplay between coherence and saturation

We solve the CCFM equation numerically in the presence of a boundary condition which effectively incorporates the non-linear dynamics. We retain the full dependence of the unintegrated gluon distribution on the coherence scale, and extract the saturation momentum. The resulting saturation scale is a function of both rapidity and the coherence momentum. In Deep Inelastic Scattering this will lead to a dependence of the saturation scale on the photon virtuality in addition to the usual x-Bjorken dependence. At asymptotic energies the interplay between the perturbative non-linear physics, and that of the QCD coherence, leads to an interesting and novel dynamics where the saturation momentum itself eventually saturates. We also investigate various implementations of the"non-Sudakov"form factor. It is shown that the non-linear dynamics leads to almost identical results for different form factors. Finally, different choices of the scale of the running coupling are analyzed and implications for the phenomenology are discussed.


Introduction
The linear QCD evolution equations at small-x, such as BFKL [1] or CCFM [2][3][4][5][6], lead to a very rapid growth of the gluon density. Once the density is too large however, one needs to take into account the non-linear effects. The BFKL evolution is modified in the high density limit by the so-called Color Glass Condensate (CGC) [7,8] (for a review see [9]) formalism where the strong color fields are treated semi-classically, and where their presence leads to the saturation of the gluon emission rate 1 . The BFKL evolution equation is then generalized to either the Balitsky-Kovchegov (BK) [10,11] equation or to the more complete JIMWLK [12,13] equation. The energy dependent momentum scale characterizing the transition from the weak to strong field regimes is called the saturation momentum and is denoted Q s (x).
What makes CCFM different from the standard BFKL-type picture of the small-x regime is the property of the so-called coherence of emissions. This is quantum mechanical in origin and comes from the interference in the emission of gluons from different partons. This leads to what is termed angular ordering which, as a fundamental feature of QCD, is present in both time-like and space-like evolution as well as at both small and large x. In CCFM the effects of angular ordering are encoded in an additional momentum scale, usually denoted p, so that the unintegrated gluon density A now depends on three variables: Bjorken x, the momentum k, and the additional momentum p. We will denote the unintegrated gluon density by A(x, k, p). In comparison, the BFKL gluon density, and similarly the density in the CGC formalism, depends only on x and k, that is A = A(x, k).
The saturation scale Q s is the momentum at which the gluons in the hadron start to overlap with each other. This leads to the increased rate for the recombination of gluons. Q s can be defined as the momentum for which A is a given constant, say of the order of 1. This is a natural definition since the gluon density can be defined as a phase space occupation number 2 of gluons per unit rapidity with given transverse momentum k. This automatically implies that Q s is x dependent. However, it is then clear that in CCFM the same definition immediately leads to an extra scale dependence of Q s , since now the requirement of constant occupation number, A(x, Q s , p) = const, leads to Q s ≡ Q s (x, p). One of our main objectives in this paper is to study this novel scale dependence of Q s by numerically solving the CCFM equation in the presence of a boundary condition enforcing non-linear corrections.
The additional momentum variable p in CCFM is determined by the kinematics of the hard scattering. More precisely, the hard scattering determines the maximal angle of the allowed emissions in accordance with angular ordering. Conventionally one uses the squared angle ξ = q 2 ⊥ /q 2 l as the angular variable, and for example in a DIS process with n gluon emissions, angular ordering implies that ξ 1 < ξ 2 < · · · < ξ n <ξ whereξ is the maximal angle determined by the hard scattering. The transverse momentum p is then defined viaξ = p 2 /(x 2 n s) where x n is the energy fraction of the gluon entering the hard scattering, and s is the cms energy. The variable p is related to Q in DIS, roughlȳ ξ ≃ Q 2 /(x 2 s). This means that the saturation scale from CCFM in DIS kinematics will depend on the photon virtuality Q in addition to the Bjorken x. In the usual picture of the small-x physics, however, Q s is an intrinsic property of the target hadron, and is independent of the probe and its virtuality. The difference now is that the gluon occupation number in the hadron depends additionally on the probe. That Q s is not a property of the target alone but also depends on the probe and its virtuality is certainly what one would generally expect from quantum mechanics. In the end Q s is determined by the scattering process, to which both the probe and the target enters, and therefore we should expect Q s to depend on both. In that sense the extra scale dependence of Q s is most natural.
Unlike BFKL, the proper non-linear generalization of CCFM is not know yet, and we shall not present it here either. In [14,15], however, a method was proposed which effectively implements saturation and unitarity corrections within a generic linear evolution equation, and the method was in particular applied to BFKL and CCFM. The idea is simply based on enforcing a boundary condition in the linear evolution which prevents the gluon phase-space occupation numbers to grow too large, and the procedure is an extension of a strategy originally introduced in connection with analytic studies of BFKL evolution in the presence of saturation effects [16,17]. The method was demonstrated to work very well for BFKL in [14] where the solutions of BFKL in the presence of the saturation boundary condition were compared to the solutions of BK, and it was shown that the method reproduces successfully the results from BK for any x and for both fixed and runningᾱ s .
The method was then applied to CCFM in [15] and Q s was extracted. However, in order to simplify the numerical procedure, only the case k ≪ p was studied in which case one can neglect the dependence on p so that Q s depends only on x. This is the region that leads to another formulation of CCFM called the Linked Dipole Chain (LDC) model [18] which has some advantages compared to CCFM such as projectile-target symmetry. In this paper, we will keep the full p dependence of the evolution, and we will show that this dependence, which as mentioned above is related to coherence in the emissions, combined with the physics of saturation gives rise to a very interesting and novel dynamics which eventually leads to the saturation of the saturation momentum Q s itself. This happens when the solution to A "stalls" in some region of k, meaning that the evolution in Y is almost completely stopped. This behavior is in contrast to the standard picture of the non-linear dynamics where Q s grows indefinitely with Y .
Generally we shall see that the shape of Q s and that of A close to the saturation regime is in CCFM rather different than in BFKL. Also some characteristic features of the standard small-x evolution such as geometric scaling [19] is modified in CCFM. On the practical level, however, it seems that the new behavior of Q s , whereby it eventually saturates as a function of Y for given p, sets in very late, especially with a running coupling, so that it seems it would be be very hard to observe at collider experiments. On the other hand it will still be important to take into account the variation of Q s with p also for phenomenological values of Y , especially when p is small.
In this study, we restrict ourselves to the version of CCFM containing only the "hard" emissions associated with the 1/z pole in the gluon splitting function. In addition to these, in CCFM there are also the "soft" emissions which are associated with the opposite 1/(1 − z) pole 3 . Within the accuracy of the formulation of CCFM however, these emissions are exactly "probability conserving", meaning that they are exactly compensated by the associated virtual corrections encoded in the so-called "Sudakov" form factor. Thus at least formally, it should make no difference for the solution of A whether one includes the soft emissions or not. The difference would be only in the exclusive final state, a correct description of which requires the soft emissions. In practice, however, such a formal statement carries little significance, and the inclusion of soft emissions can have substantial effects on A. Thus it is our aim to also include the soft emissions in our procedure. This is, however, somewhat challenging numerically due to the fact that the real and virtual terms are separately divergent, and so one must introduce a cut-off which complicates the numerical procedure. Therefore the inclusion of the potentially important soft emissions and the associated Sudakov form factor is postponed to a future work.
The 1/z emissions are accompanied by virtual corrections which are encoded in the so-called "non-Sudakov" form factor. One can in the literature find different expressions for this form factor. In this paper we shall study both of the common expressions for the form factor. We shall see that, generally, the different choices of the non-Sudakov form factors lead to rather different results in the linear evolution, which is due to the reason that the differences which occur in the softer region are enhanced unless there is some mechanism to cut off the growth there. As we will demonstrate, once non-linear evolution is included via the saturation boundary condition, then the differences are to a large degree removed. The small differences remaining when Y is extended to very large values can can be removed by scale factors independent of Y and p. Thus saturation introduces a certain universality in the evolution as was already observed in [15].
Related to this, it is important to check the sensitivity of the results to different prescriptions of the boundary condition. Generally, one can consider different boundary conditions or, for a given type of boundary, different values of the boundary parameters. In [15] it was reported that the different choices lead to different normalizations of the solutions, without altering the shape or the Y dependence. Since we here keep the p dependence, the question is more subtle. We will show that the differences are small and that they can generally be removed by Y and p independent scale factors. Some differences can remain in the region where Q s starts to saturate (the region where A "stalls"), but what we find is that this behavior itself, i.e the saturation of Q s , is independent of the precise application of the boundary. What changes is simply at what Y , for a given p, this behavior sets in. To know the exact details of the behavior in the saturated regime one would need to fully formulate the proper non-linear generalization of CCFM which is beyond the scope of this work. It seems, however, rather reasonable to expect that this interesting dynamics where, as a result of the combination of saturation and coherence, Q s eventually saturates will remain true also once the proper non-linear generalization of CCFM is known.
The structure of the paper is the following: in the next section we present a brief overview of the CCFM formalism and show the numerical results for the solution in the case when the full dependence on the scale p is included. We also discuss different versions of the running coupling prescriptions. In section 3 we provide semi-analytical results for the case of CCFM with the boundary condition. In particular, we discuss the dependence of the saturation scale on the rapidity and p. In section 4 we briefly describe the method of implementation of the saturation boundary. In section 5 we present the complete numerical analysis of CCFM in the presence of the saturation boundary. We analyze different prescriptions for the non-Sudakov form factors and the running coupling, and we extract the rapidity and p dependence of the saturation scale. Finally, in section 6 we state our conclusions and discuss the relevance of the results for the phenomenology.

Overview of CCFM
Our aim in this section is to provide a brief overview of CCFM, focusing on the points which shall be relevant for our analysis. It is possible to find different versions of the non-Sudakov form factor in the literature but the origin of this form factor is hardly ever discussed. It will therefore be important to recall the derivation of this form factor and understand the motivation behind the formulas. Moreover, we shall in addition to the standard choice of the scale of the running coupling also consider a different choice. The most comprehensive overview is to be found in the original papers in [2][3][4][5][6], and a simplified but rather detailed discussion on CCFM is also offered in [15]. The gluon ladder generated by the CCFM evolution is represented in figure 1. In CCFM the successive emissions along the ladder are ordered in their angle. As mentioned in the introduction, it is customary to define ξ = q 2 ⊥ /q 2 l as the angular variable (with the direction of the incoming proton defining the longitudinal direction), and angular ordering implies that ξ i > ξ i−1 . Additionally there is a maximum angleξ, set by the kinematics of the hard scattering, which limits all emissions.
In CCFM, the emissions are divided into hard and soft emissions, and the importance of this classification is that it offers significant simplifications in the formulation. The precise definition is that the hard emissions are those emissions ordered in both energy and angle, while for each soft emission there exists another emission with bigger angle and larger energy fraction. It follows that, the hard emissions are those associated with the 1/z pole in the gluon splitting function while the soft emissions are associated with the opposite 1/(1 − z) pole. Following the above definition one can associate with each hard emission k a cluster of soft emissions S k which satisfy ξ k−1 < ξ < ξ k and y < y k where ξ k−1 is the angle of the previous hard emission. The label "soft" has a two-fold meaning. Firstly the fact that each of these emissions takes away very little energy, because of the 1/(1 − z) pole, and secondly that all soft emissions in the set S k satisfy q < q k . The classification of the emissions as being hard or soft is also related to another classification of the ladder in terms of so-called k ⊥ -changing and k ⊥ -conserving emissions respectively. What is meant by this is that there are certain real emissions which leave the virtual propagator momenta unchanged in an approximate sense. This latter classification is crucial for deriving the 'standard' version of CCFM which is used in all practical applications.

Integral equations and the virtual form factors
The CCFM integral equation, including both the soft and the hard emissions, is given by (2.1) The momentum variable p is defined viaξ = p 2 /(x 2 n s), and k ′ = |k k k + (1− z)q q q|. The momentumq is the rescaled momentum of the real gluon, and is related to the true momentum q byq = q/(1 − z). In line with all the approximations made in order to derive this equation one can for the 1/z pole setq = q (which is what we have already done in ∆ ns (z, k, q)). Here ∆ s is the Sudakov form factor which screens the singularity of the 1/(1 − z) pole, while ∆ ns is the "non-Sudakov" form factor which depends also on the momenta k and it corresponds in BFKL to the gluon Regge factor. The rescaled coupling constantᾱ s is defined as α s N c /π.
In figure 2 we represent pictorially the components of equation (2.1). Here the horizontal axis represents the inverse energy fraction of the gluons, the further to the right a gluon is located the less energy it has, while the vertical axis represents the transverse momenta, the higher up a emission is located the larger its transverse momentum. The maximal angleξ is represented by drawing a line of constant angle from the point (x n , p); this is the diagonal line going through the black dot representing p (p itself is not the momentum of any gluon (real or virtual) but we still denote its location by a dot in the figure). In this particular case, k < p. The "Sudakov" form factor ∆ s (p, zq) is then given by exp(−ᾱ s S) where S is the phase space area which is the diagonally shaded region in the figure. This is the region where further emissions would have been possible. Similarly the "non-Sudakov" form factor ∆ ns is given by exp(−ᾱ s A) where A is the area of the vertically shaded region. In the figure we also show the next backward step in the evolution which is obtained from (2.1) when A in the RHS is expanded iteratively. As we can see from the figure, the region S contributing to ∆ s diverges when extended down indefinitely. Thus this divergence is cut by a horizontal line in the figure, which means a soft cut-off in the transverse momentum. Notice that this divergence is the one arising from the 1/(1 − z) pole because as z → 1, the gluon carries infinitely small energy fraction, and in the figure it would then be located precisely in the region where the diagonally shaded region is extended down indefinitely (the fact that this region is bounded by the diagonal lines is due to angular ordering).
The soft emissions are such that they are exactly compensated by the factors ∆ s . This means that the regions over which the Sudakov form factors receive contribution are precisely the regions to which the soft emissions are confined, i.e. the clusters S k mentioned above. Summation over all possible soft emissions in each S k leads to a factor exp(+ᾱ s S k ). Each such factor then multiplies to unity with the corresponding Sudakov form factor exp(−ᾱ s S k ). One is consequently left with the hard emissions only, and the integral equation reads This equation is thus formally equivalent to (2.1). In practice, however, one can expect to find different solutions. The numerical implementation of the soft emissions faces the difficulty that one must cut the z → 1 singularity by a momentum cut-off as in figure 2. This introduces a numerical uncertainty which complicates the procedure. Therefore, even though their inclusion is important and interesting, we will postpone their treatment to a future work and in the present concentrate only on equation (2.2). As mentioned above, the non-Sudakov form factor ∆ ns can be written as exp(−ᾱ s A), where in figure 2, A is the vertically dashed region bounded by k and the angle of q, and the energies of q and k. In CCFM, the form factors ∆ ns and ∆ s are derived out of the socalled eikonal, S eik , and non-eikonal, S ne , form factors which sum up to all orders the virtual corrections associated with the eikonal and non-eikonal emission vertices respectively [2][3][4][5][6]. As shown in [15], for each emission k the product of the associated form factors S eik (k) and S ne (k) equals the BFKL form factor ∆ BF KL (k), that is S eik (k) · S ne (k) = ∆ BF KL (k). The definition of ∆ ns for a given emission q is that it is the product between S ne and a portion of the eikonal form factor which covers a region bounded from below by the angle of q and from above by the maximal angleξ (this is the region left over after multiplication with a soft factor coming from the resummed real emissions). In figure 3 this is the sum of the shaded regions B 1 and C 1 . The non-eikonal form factor S ne itself is such that it covers a region bounded from below by k while from above it is again bounded byξ, this is region B 1 in the same figure. What is important however is that the sign in the exponential factor is reversed, so that it is positive. 4 The non-Sudakov is then given by In this case, however, we assumed that k ≥ q, see figure 3. Now, if k < q we can have two different situations, depending on whether or not k is entirely below the diagonal line through q in the figure, that is whether k < zq or k ≥ zq. In these both cases we define a new region A i shown in figure 4. This is a region bounded from below by k and from above by the angle of q. When k > q as in figure 3 of course A 1 = 0.
On the left plot of figure 4 we show the case when k ≥ zq, and the definition of ∆ ns then yields  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000  00000   11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  11111  If A 2 is bigger than C 2 then we see that ∆ ns is larger than unity. In case k < zq we instead have (in this case C 3 = 0) Thus ∆ ns exceeds unity when A i ≥ C i , and this actually happens when k 2 < zq 2 . All the results above give Therefore if the kinematical constraint k 2 ≥ zq 2 is assumed to be valid, ∆ ns is guaranteed to be bounded by unity 5 . A different version of ∆ ns was presented in [21] in which case the regions A i are simply neglected. Then one always has where Since A i is neglected, we are thus guaranteed that ∆ ns ≤ 1. Thus in a sense this has a similar effect to assuming the kinematical constraint to hold, but one should still keep in mind that the effect is not the same because one has no apparent reason to simply neglect the regions A i (this amounts to neglecting soft emissions above k, see also discussion below). What the kinematical constraint ensures is that C i ≥ A i . In this paper, we will not include the kinematical constraint (it was included in [15]), which must be included in the real emissions as well, but we shall rather study the differences between using the two different non-Sudakov form factors, (2.6) and (2.7).

Scale of the running coupling
Above we assumed a fixed coupling which is consistent with formulation of CCFM which was only performed at LO where the coupling is fixed. Although the effects of the running coupling are formally of next to leading order, it is nevertheless for any phenomenological application extremely important to take into account the running of the coupling. The standard procedure is to let the coupling run with k in the 1/z pole and in ∆ ns which is similar to what is usually done in LO BFKL. When both hard and soft emissions are included then q is usually chosen as the relevant scale forᾱ s in the 1/(1 − z) pole and in ∆ s (this is for example done in CASCADE as well as in the SMALLX program [22]). Notice that ∆ s multiplies also the 1/z pole so with this choice one is for the hard emissions using q as scale for the Sudakov form factor while k is used for the non-Sudakov and the real emissions.
As k is the larger scale in ∆ ns it may seem natural that it is chosen as the scale for the running coupling. However, once a scale is chosen for a certain class of emissions, say the hard emissions, then that scale should be consistently applied everywhere because otherwise some cancelations which are important for the consistency of the formalism may be violated. Also, even though k seems to be the larger scale in ∆ ns we have seen above that actually in the derivation of ∆ ns , k first enters as the lower scale in S ne .
In [22] it is argued that the choice to have q and k as scales simultaneously can be derived from having the single choice (1 − z) |k 2 | (where k 2 is now the four momentum squared), since this choice reduces to k and q when z → 0 and z → 1 respectively. One should, however, again remember that both ∆ ns and ∆ s are derived out of S eik and S ne , and while the four momentum of k could be argued as scale in the latter, the former resums the graphs where virtual gluons are emitted and absorbed in between the real gluons so it is hard to see why the scale of k should also be used here. We find it more plausible that instead the momenta q of the virtual gluons are used in S eik . But if this choice is made, then to derive ∆ ns and ∆ s in the first place, one has to choose q as scale in S ne too, since otherwise necessary cancelations between S eik and S ne would not work. Consequently this leads to using q in both ∆ ns and ∆ s , and therefore also in all the real emissions.
If one uses k as scale in ∆ ns then obviously the formula in (2.6), and similarly the corresponding one in (2.7) is the same as in the fixed coupling case. If on the other hand one uses the virtual momentum q ′ as scale forᾱ s , then the formula is changed. We use q ′ as scale in (2.7) and derive a new formula which is however somewhat lengthy so we do not present it explicitly. We shall use the one loop result for the running coupling, frozen below some scale k 0 , Thus for the running coupling case we shall consider the following alternatives: One in which we use k as scale and ∆ ns given in (2.6), one in which we use k as scale and ∆ ns in (2.7), and the last case where q is chosen as scale inᾱ s and the corresponding modified ∆ ns is used. We shall note here that, in the BFKL equation the coupling also starts to run at the next-to-leading logarithmic level (in ln 1/x) only [23][24][25][26][27][28]. The different choices of scales formally lead to the differences at next-to-next-to leading level. They are however not negligible numerically. The form of the NLL correction suggests that the most natural scale is that of the transverse momentum squared of the emitted gluon, see for example [29,30]. This confirms the choice of (2.9) as a more natural choice of scale. We will present here numerical results for the solution to the linear CCFM evolution. We start with the comparison of the two versions of the non-Sudakov form factor (2.6) and (2.7) in the case when the coupling is fixed, and set to be equal toᾱ s = 0.2. In figure 5 we present the value of the unintegrated gluon distribution A(x, k, p) obtained from CCFM as a function of the transverse momentum k, and for the fixed value of the momentum p = 10 GeV (left plot) and p = 200 GeV (right plot). Different curves in increasing order correspond to increasing values of the rapidity. The initial condition for the equation was set to be with µ = 1 GeV. It is clear that the 'non-regular' form-factor given by (2.6) gives a significant contribution in the low momentum regime due to the fact that k 2 < zq 2 so that the expression in the exponent changes sign. This gives a substantial enhancement at low values of k. On the other hand it is also clear that both prescriptions for the formfactors coincide in the large k regime. One observes that the slope of the k distribution becomes steeper when k becomes larger than p, meaning a suppression of the momenta due to angular ordering. Momenta larger than p are nevertheless allowed. This change in the slope in k is crucial for the analysis with the boundary and the behavior of the saturation scale as we shall see later.

Fixed vs running coupling
In figure 6 we perform the comparison between the running and fixed coupling scenarios. The running coupling is regularized at k 2 0 = 0.7 GeV 2 with frozen scenario. As expected the running coupling gives the higher contribution in the region of the infrared momenta. It falls then below the fixed coupling at high momenta but such that k < p, as the running coupling value becomes smaller than the fixed valueᾱ s = 0.2. The effect is most pronounced for the calculation with p = 200 GeV. We observe though that in the region where k > p, the scenario with the running coupling gives results which tend to be larger. One needs to take into account the fact that in this region we have stronger relative suppression from the   non-Sudakov form-factor (we shall discuss this more below). The value of the form-factor is smaller for the larger value of α s , which is for the case of the fixed coupling. This is why the running coupling calculation tends to be less suppressed in the asymptotic regime of the large transverse momenta.

Comparison between different versions of the running coupling
Next let us turn to the implementation of different scales in the running coupling. In figure  7 we show the comparison between the standard results and those obtained using q as scale. The differences between two choices are in general not too large, with the choice of q giving smaller values for the solution. In the case when p ≫ k there are more and more chains with q ≫ k contributing, and so typicallyᾱ s (k) >ᾱ s (q). Therefore the solution with q as a choice of the scale is slightly lower.
On the other hand, when p is small (∼ 1 GeV) we see that the differences are larger, and especially the Y dependence and the slope of A seem to be modified. In this regime the phase space is strongly constrained by the maximum angle allowed for the emissions. The phase space restriction applies to the real momenta q, while the virtual momenta k are not directly restricted. Then in an event where k ≫ p we would expect that generally k ≫ q as well (except early in the chain). This implies thatᾱ s (k) <ᾱ s (q) and one could therefore expect to see a faster growth when q is chosen as the scale. However, the opposite behavior is observed. This is because a larger coupling also implies that there is a larger suppression from ∆ ns . We saw a hint of this behavior already in figure 5. When k > p the phase space for the real emissions is more strongly constrained, on the other hand (2.6) or (2.7) are independent of p, so even if the real phase space is strongly constrained, the form-factor is not. This implies the possibility that a larger coupling in this regime actually slows down the evolution due to the large suppression from ∆ ns .
To further demonstrate that this is indeed the case we show in figure 8   For clarity of demonstration we show only solutions for Y = 8, 12, 16. Of course in this case sinceᾱ s = 0.3 >ᾱ s = 0.15 for any k, in the region where k p the solution with the larger coupling indeed grows faster (the phase space restriction is not so important here). However, we then clearly see that as k > p, the suppression from ∆ ns , which is much larger whenᾱ s = 0.3, starts to dominate and therefore the solution with the larger coupling is strongly suppressed. Consequently in the region k ≫ p, the solution with the larger coupling actually grows slower.

CCFM in the presence of saturation: analytical results
Having reviewed some of the basic concepts in CCFM we now move on to the question of implementing saturation effects and finding the relevant saturation scale Q s . The implementation of saturation effects will be done using the same approach as in [14,15], but most importantly because in this work we do keep the full dependence of A on the additional scale p, we shall discover many interesting new results not noticed before. The main implementation is numerical and all the numerical results are to be presented in section 5. In this section we will analytically study the general behavior of Q s .
To study the saturation scale it is convenient to take the Mellin transform over the momentum k. Following [31] we shall write the Mellin representation of A as where Y = ln 1/x, and we define ρ = ln(k 2 /Q 2 0 ) with Q 0 an arbitrary scale. We will only consider the fixed coupling case in the analytic formulas. The function H encodes the dependence of A on the parameter p. What we know generally is that H ≈ 1 when p/k ≫ 1, since the restriction coming fromξ becomes irrelevant in that case. This, however, does not imply that angular ordering is not important anymore. More precisely it means that the ordering in the final step of the evolution is automatic when p ≫ k, but in the previous steps of the evolution it is still important. This will reflect itself in the characteristic function ω which in the limit p ≫ k is given by To find the function H one can differentiate the integral equation (2.2) with respect to p as done in [31]. One then finds that where as before k ′ = |k k k + q q q|. Assuming that Y is very large and the solution is in the asymptotic region one can neglect the first theta function, which is what is done in [31]. However, one should then notice that with the non-Sudakov in (2.6) which is allowed to exceed unity, one would run into a divergence. This is so because when Y → ∞, even if p is finite, the phase space for the emissions contributing to the scattering becomes infinitely large, and there is a UV divergence appearing when q → ∞. On the other hand such a divergence does not appear if the restricted ∆ ns in (2.7) is used instead. The reason for this is that the choice in (2.7) effectively implies that k plays the role of a UV cut (but only for the soft emissions). If the restriction fromξ is removed (as when Y → ∞) then the form factors S eik and S ne need to be regularized by some UV cut-off. When these factors are multiplied together, because of the different signs, they cancel each other in the region above k, and thus the dependence on the UV cut-off vanishes. However, the real emissions need still be regularized by the UV cut-off, since otherwise the inclusive summation over the soft emissions for the given hard emission q would give a divergence as q → ∞. The fact that the regions A i in figures 3 and 4 are neglected in the derivation of (2.7) means that one is cutting off all real soft emissions above k which is why k plays the role of the UV regulator 6 for the soft emissions. Thus for simplicity we shall use this form factor in the following, though a proper renormalization procedure for all the emissions would probably be called for. After these remarks, let us now return to the equation for H above. The saturation of the solution can be analyzed in the two different regimes of large and small p.
The case when p > k: In case p > k, since q > p (because of the theta function in (3.3)), we also have that q > k. We will then approximate k ′ ≈ q. Equation (3.3) is then easily solved to give where in the second equality we neglected the term exponential in Y , because it is very small when Y is large. We have also defined ρ p ≡ ln(p 2 /Q 2 0 ). Consistency of the solution requires that We see that when (ρ p − ρ) ≫ 1, that is when p is very large, H reduces to unity. Since (ρ p − ρ) is large we can rewrite H as The Mellin representation then reads .

(3.7)
To find the saturation scale, one can essentially follow a line of constant density [16,17] which means that the relevant saddle point for the saturation problem is determined by the equations (defining ρ s ≡ ln(Q 2 Here all derivatives are with respect to γ and the second equation is the usual saddle point condition while the first equation states that the integrand around the saturation scale is constant and of order 1. On the other hand we recall that in the saturation problem for BFKL the corresponding saddle point equations read What we notice, besides the fact that the equations in CCFM look much more complicated than in BFKL, is that equations (3.8) and (3.9) imply that the saturation anomalous dimension γ s , and thus also the characteristic function ω s , depend generally on p and Y . In BFKL on the other hand, γ s is a pure number and ω s is independent of Y which we see from the corresponding equations. As ρ p − ρ s → ∞, we see that the structure of the CCFM saddle point equations reduce to the BFKL ones, though at finiteᾱ s , γ s and ω s are still different than in BFKL. Whenᾱ s → 0, the CCFM solution reduces completely to the BFKL one.
From equation (3.8) we can find ρ s . Let ρ (0) s be the solution when ρ p − ρ s → ∞, that is Since the general form of ω which depends on p is complicated, we will in the following simply make the approximations to set γ s = γ Thus we find that which is valid as long as ρ p − ρ (0) s ≫ 1. The saturation saddle point is determined by as in BFKL. Of course this is the lowest order approximation (lowest order with respect to the corrections in p, but it is valid at finiteᾱ s ) but since the more general case is rather complicated we shall use this formula. In the later section we will determine the full saturation momentum via the numerical solution.
In order to be able to use the formula (3.15) one needs to have a solution for the characteristic function in CCFM. This is essential for obtaining ω  s . Unfortunately it is only possible to extract the characteristic function numerically, as was done in [31]. To It is evident that for small rapidities the growth of the saturation scale is dominated by the term exponential in rapidity and the characteristics is similar to the BK equation (or BFKL with the saturation boundary). However for fixed value of scale p the dependence of ρ s on rapidity becomes nonlinear. The smaller the value of p is, the lower is the critical value of rapidity for which the behavior starts to become nonlinear, and there is an additional dependence on p. The formula cannot be trusted in the region where ρ s falls off with rapidity. This region is outside the validity of the original approximation in which equation (3.15) was derived. One would need to compute higher order corrections and possibly resum them in order to predict the asymptotic behavior of Q s with Y .
Thus the leading term for the saturation momentum when p is very large has the same structure as in BFKL, with the "speed", −ω ′ s , of ρ s being much slower, around 0.65 for CCFM as opposed to around 0.98 for BFKL at LO andᾱ s = 0.2. The corrections coming from finite p slow down the saturation scale and from (3.15) we see that the p dependence in Q 2 s enters as a power inside an exponential. Summarizing, we observe that the corrections from the phase space restriction in CCFM introduce an increasingly complicated shape of Q s depending both on p and Y when going beyond the BFKL regime, i.e. p ≫ k.
Using the results above one can also study the behavior of A near saturation, that is when 0 < ρ − ρ s 1. This is done by expanding the function in the exponential of the Mellin representation to O((ρ − ρ s ) 2 ). In BFKL this leads to The BFKL result features the well known characteristics of the small-x evolution. One is the property of geometric scaling which is the fact that the Y and k ⊥ dependence enters the gluon distribution in the combination k ⊥ /Q s , that is via ρ − ρ s . Indeed as long as ρ − ρ s ≪ 2ω ′′ s Y the dominant dependence of A is on ρ − ρ s . A more careful treatment performed in [16] modifies this solution so that the Y dependent pre-factor coming from the Gaussian integration is replaced by a factor (ρ − ρ s + ∆) for some number ∆. In this case geometric scaling is even more obvious. The second property is the diffusion in k ⊥ space. This is governed by the Gaussian and the diffusion radius grows as 2ω ′′ s Y . One can perform the exact same manipulations also in (3.7). However the resulting expression is rather long and not very transparent so we do not write it down here. Again in case ρ p − ρ ≫ 1 the structure reduces to that in (3.18). If the p dependence is kept, both of the two characteristics mentioned above are modified. For example, diffusion is reduced due to angular ordering, and the diffusion radius obtains a complicated dependence on ρ p − ρ. Geometric scaling is also not transparent anymore, because of the additional p and Y dependence. Besides, generally we would need to take into account the p and Y dependence of γ s and ω s which further breaks the scaling behavior. One may find a modified scaling behavior even if it is not as transparent as in the BFKL case, but given the complexity of the solutions we shall refrain from pursuing this question analytically.

The case when p < k:
The opposite region where p ≪ k is more complicated, but also turns out to be much more interesting. In this case we can no longer use expression (3.2) for ω, and the p dependence is stronger. Note that in the argument above, which follows the one in [17], saturation is not really explicit. It is understood that the solution represented by the Mellin transform is valid only as long as one is above Q s , but the effects of saturation are kept implicit. Thus one could actually get the same result if one simply took linear BFKL and followed a line of constant amplitude, though of course in the linear case this is rather ad hoc since there is nothing special about any particular line of constant amplitude, nevertheless the answer for Q s would be the same. Therefore the leading asymptotic Y behavior of Q s extracted in this way from BFKL, without actually explicitly applying saturation, is the same as that extracted from the non-linear BK equation. To obtain correctly the non-leading Y behavior, however, one has to go beyond this approach but what was shown in [16] is that one can get the correct answer from the linear equation if the effects of saturation are introduced via the boundary prescription.
In CCFM the effects of saturation will be seen to be rather dramatic. From the numerical analysis we shall see that the phase space restriction coming from p, which has its origin in the coherence of the QCD radiation, when combined with saturation, which imposes a restriction in the region below Q s , implies that the saturation scale itself saturates at a region Q s ≫ p. This behavior is missed by the naive analysis above since one must explicitly apply saturation to possibly see it. A treatment like in [16] may lead to an analytic understanding of this novel behavior, but given the complexity of the equations involved this may be very difficult to pursue. We shall therefore not apply the analytic treatment to this case, even though finding H itself is easy, as done in [31]. Instead we move on to the numerical treatment.

Description of the boundary method
The boundary method applied to the linear BFKL evolution, was shown in [14,15] to successfully reproduce the results obtained from the non-linear BK equation. The saturation momentum Q s and the shape of the gluon distribution for k ≥ Q s were correctly reproduced by this method. The basic idea is to prevent the gluon occupation number from growing too large. First, we define a line of constant occupation number via the condition where the number c is of order 1 or smaller (but not much smaller). As obvious from the definition, ρ c depends on both Y and p. In the linear evolution, A will continue to grow indefinitely beyond ρ c when increasing Y . We then enforce a boundary condition preventing such an unrestricted growth beyond the critical line ρ c . In practice this is done by applying the boundary for all momenta ρ ≤ ρ c (Y, p) − ∆. Here ∆ is a given number.
Strictly speaking in CCFM ∆ should not be a pure number but should depend on Y and p. However, given that explicit analytic solutions of CCFM are not available, it is not really feasible to analytically find the possible dependence of ∆ on Y and p. We will therefore for simplicity let ∆ be a constant as in BFKL. The parameters c and ∆ can be viewed as free parameters although in BFKL it is known that ∆ ∼ ln(1/c). In particular the results should not depend on the precise values of these parameters, apart from differences in normalization. Our default choice will be to set ∆ = 2 and c = 0.4. In [14,15] the default choice of the boundary condition was to let In this case one has a totally absorptive boundary, and this was indeed also the choice in the analytical analysis in [16]. This condition is convenient for the analytical analysis in the BFKL case because then the problem can be formulated as diffusion in the presence of an absorbing wall [16]. The precise choice of the boundary condition, should (at least asymptotically) not make any difference beyond a normalization factor. Indeed , it was explicitly demonstrated in [14,15] using the numerical analysis that the method is robust for any Y (thus not only asymptotic ones), so that different choices still give similar results. In this paper we shall not be using the above condition, but we will rather investigate different alternatives.
One possible choice of boundary is to let A grow logarithmically beyond the lineρ c . Then we let This is similar to the behavior observed in the BK equation.
Another choice is to simply freeze the gluon distribution at the critical point, Once the boundary has been implemented, one can extract the saturation momentum from the evolution. Q s can be defined as a line of constant density. One can define ρ s as any line of constant density between ρ c andρ c . The precise choice affects only the normalization of Q s .
We will investigate both choices of the boundary condition. This is because one does not really know the details of the non-linear dynamics in the CCFM, and we want to check the robustness of the obtained results with respect to the change of the boundary condition.

Fixed coupling results
In figure 10 we show the solution to the CCFM with the saturation boundary implemented in the form of freezing as in (4.4). Here, the coupling is fixedᾱ s = 0.2. We compare solutions obtained using (2.6) and (2.7) respectively. We see that once the saturation effects are included, the differences between the solutions using the different form-factors are not that large, as the sensitive infrared region is regularized by the saturation boundary. The dependence on p can be examined by comparing the left and right plot in figure 10. We observe that, similarly to the linear case, the distributions change slopes when k become comparable to the external scale p. For k > p, the gluon distribution is effectively cut off by the scale p.
In order to analyze the asymptotic behavior of the solutions we extend the rapidity range up to Y = 50. This is shown in figure 11, where we use the non-Sudakov form-factor from equation (2.6). We observe the interesting effect that the front of the solution no longer travels to higher momenta with uniform 'speed' governed by the saturation scale, but that it rather stalls at high rapidities. The origin of this behavior is rooted in the presence of the additional scale p. As the rapidity grows, the solution diffuses to the higher values of the transverse momenta k. At the same time the saturated region grows also with energy. The unsaturated region is eventually pushed to the region where k ≫ p, which is then cut off by the angular ordering condition. Therefore the phase space for new emissions is squeezed, due to the presence of two scales Q s and p. We have done the analysis with two different boundaries (4.3) (logarithmic growth at low momenta) and (4.4) (freezing) and it is clear that the front stops to diffuse to higher momenta in both cases. The stalling of the front therefore is a robust feature of CCFM with saturation as it does not depend on the way the boundary is implemented. We note that the complete stalling of the front (in the case presented in figure 11) occurs for momenta k ∼ 10 2 GeV, which is order of magnitude larger than the p ∼ 10. In consequence, the saturation scale Q s starts to behave differently with rapidity in this regime. In particular the value of ln Q s will no longer scale linearly with rapidity but slower than that, and eventually the saturation scale will 'saturate' itself at large Y . The complete analysis of the behavior of the saturation scale is presented in the later section.

Running coupling results
Let us now turn the case of a running coupling. In figure 12 we show the different results obtained using the two different form factors, (2.6) and (2.7). Here, the boundary condition (4.4) is used, and k is used as scale inᾱ s . As in the fixed coupling case we see that the differences are much smaller than in the linear evolution, and almost identical results are obtained. In particular we have checked that the small difference between the results stays the same up to very large values of Y (far beyond phenomenologically relevant values). We find that, apart from the region of smallest Y in the figures, the difference is a normalization factor independent of p and Y . Consequently the different form factors lead to almost identical saturation scales which we shall look more closely at in the next section. Since in the running coupling case the evolution is generally slowed down, one needs to go to higher values of Y to see the stalling of A. For larger p in fact, the Y values required for the stalling are so extremely large that we have not reached them in our numerical solution. Starting from smaller p values, however, we clearly see how the evolution is eventually slowed down to a point where A starts to stall in a region betweenρ c and ρ c (i.e. around ρ s ), and one can also see how this region then grows with Y . Increasing p we see that one needs ever higher Y to see this behavior (as in the fixed coupling case).
In figure 13 we show the extension of the solution shown in the left plot of figure 12 (that is when p = 1 GeV) to asymptotic values of Y between 46 and 70 units. In this case only the solution using (2.7) is presented, the solution using (2.6) looks essentially the same. Compared to the fixed coupling case, the region where A stalls is much smaller, and it also grows slower. Therefore the saturation of the saturation scale will set in much later (to be shown in the next section). To be able to exactly determine at what Y the stalling of A starts to become visible, we would need to know the full non-linear formulation of CCFM. The application of different boundary conditions seem to, however, indicate that this behavior sets in rather late, as manifest in figure 13. Moreover there are additional effects not taken into account here (such as energy conservation) which will further slow down the evolution.
Let us now turn to the implementation of different scales inᾱ s . In the linear case we observed the fact that generally the choice of q gives a slower evolution compared to the default choice k. In figure 14 we show the comparison between the results obtained from the two different scales. For the larger p values the differences between the two choices of scales are very small. For smaller values of p the differences are more pronounced, the solution with scale q being consistently smaller than with the choice of scale k. As explained earlier in section 2.3.3, this effect is due to the stronger suppression from ∆ ns when q is selected as scale.
This, however, brings us to another point. Namely the fact that one should actually not allow any contribution to the non-Sudakov form factor from the region which lies outside the region set by the maximal angleξ, or equivalently by p. To understand this, remember the derivation of ∆ ns shown in figures 3 and 4. The regions B i in that case were bounded by the maximal angleξ. When k ≫ p, it can be that region B i partly, or completely, lies outside the region bounded byξ. In that case also parts of region C i will be located outside the allowed region (however, C i cannot lie completely outside the allowed region). In reality, however, one should not include any contribution from those parts of C i and B i (actually since k > q only figure 3 and the corresponding regions B 1 and C 1 are relevant). Thus one should in equation (2.3) actually only include the part of region C 1 which is bounded byξ. By not doing this one is generally overestimating the suppression from ∆ ns when k ≫ p. The non-Sudakov form factor will then be modified, and it will depend not only on k, q and z, but also on p and the previous values of z i in the chain. Unfortunately this also implies that it is no longer possible to write an integral equation as in (2.1) or (2.2) because one needs to keep track of all z i values in each ∆ ns during the evolution (this is because ∆ ns then depends not only the vertex k ′ → k + q but also on earlier vertices). One can only write down the gluon distribution in an unfolded form, explicitly summing it over all possible chains.
Thus when using the scale q instead of k we do find some differences in the evolution. On the other hand, for reasonable values of Y these differences are rather small, and for p values around 10 GeV and higher they are negligible. We will also see in the next section that the respective saturation momenta are quite similar, despite the noted differences.
We have also noted that when k > p one should actually modify ∆ ns to be consistent. However, the implementation of this modification is probably in practice only feasible in a MC simulation where one can keep track of the history of the entire gluon chain.

Fixed coupling results
From the solutions to the gluon distribution one can extract the saturation scale Q s . We have already commented on the definition of Q s in section 4. The exact value of the constant defining Q s only affects its normalization which we cannot control anyhow, the p and Y dependences are independent of the value of this constant (of course the constant should not be varied over many orders of magnitude).
In figure 15 we show the results for the two different boundary conditions (4.3) and (4.4) when the coupling is fixed and the non-Sudakov in (2.7) is used. In both cases we find a behavior consistent with the expectation in (3.15) which holds as long as Q s ≤ p. As Q s grows above p, as is most visible in the case when p = 1 GeV in figure 15, we see that the growth of ρ s = ln Q 2 s /Q 2 0 on rapidity becomes non-linear. The Y range shown in the figure is similar to what is accessed at the LHC. The absolute values of Q s should not be considered as realistic estimates of what Q s would be at the LHC (a fixed coupling calculation is not so relevant phenomenologically anyway). The two different boundary conditions give rise to slightly different Y dependences for Q s for each given p, but the differences are rather small.
The rather dramatic change in the Y dependence of Q s is visible in figure 16 where we extend the results from figure 15 to much larger Y (far beyond experimentally accessible values). As we see, the general shape of Q s does not depend on the precise boundary being used. We see that once Q s > p, the dependence of ρ s on Y becomes non-linear, and eventually it saturates to some value depending on p. Of course we cannot say exactly whether there exist a given constant C such that Q s ≤ C as Y → ∞ or whether Q s keeps growing, albeit extremely slowly. Nevertheless, it is clear that, the growth of Q s on rapidity saturates and that the dependence of ρ s on Y is highly non-linear. Interestingly, we observe that there exists a certain scaling in the sense that once Q s has saturated, the ratio between the values for two given values of p is independent of p. In other words, we find that the ratio Q s (Y, p)/p is a constant independent of p, once Q s ≫ p. The value of this constant depends on the precise boundary condition used, for example for the boundary in (4.3) we find Q s (Y, p)/p ≈ 90 when Y is very large, while for the boundary in (4.4) we instead get Q s (Y, p)/p ≈ 40.
Next we investigate the results obtained using the non-Sudakov in (2.6). In this case since the normalization of A is larger, the normalization of Q s will also be larger. We only show the results using the boundary condition (4.3), the results for (4.4) are similar. The plots for Q s are shown in figure 17 and we find exactly the same type of behavior as above. This is further seen from the right figure where the comparison to the above results are shown. As can be seen, the energy dependence of Q s is the same regardless of the non-Sudakov chosen, and we find the difference in normalization to be a factor roughly 1.35, regardless of the value of p.

Running coupling results
It is well known that the growth of Q s is slowed down by the running of the coupling. Physically it is of course consequence of the fact that the evolution leads to the diffusion towards higher transverse momenta where the coupling is much smaller. In BFKL, the linear dependence ρ s ≃ λY at fixed coupling is for a running coupling replaced by In CCFM we again expect to find a BFKLlike behavior as long as Q s ≪ p. As a consequence of the much slower growth, one needs thus in the running coupling case much higher Y values in order to reach the point where Q s > p. This implies that the stalling of the growth, and the consequent saturation of Q s is delayed to higher energies.
We start by again comparing the effects of different boundary conditions. In addition to the two boundary conditions presented above, we also try an additional boundary where we simply fix A to be a given constant behind the saturation front. In particular we choose A(Y, ρ, p) = 1, for all ρ ≤ρ c (Y, p) . (5.1) We will with this choice further demonstrate the robustness of the method.
In figure 18 we compare different boundary conditions for both of the non-Sudakov form factors. We see that the differences are minor and that they can be removed by p and Y independent scale factors. For both cases we find a difference in scaling of around a factor 1.20. We also see that Q s as in BFKL grows much slower and also that, as a consequence, the dependence on p is weaker. As in the fixed coupling case, as long as Q s ≪ p we find a BFKL-like behavior, thus in the phenomenologically relevant Y interval shown in the figure we find that for p = 200 GeV, ρ s grows proportional to √ Y . As Q s becomes comparable to p, however, this dependence is altered and becomes more complicated. Continuing to larger Y , one again finds that Q s eventually saturates. From figure 19 we also see that the difference between the results for two non-Sudakov form factors are rather minor, and turns out to be a pure scaling, around a factor 1.10, independent of p and Y . These results show that seemingly very different linear evolutions obtained from the different non-Sudakov form factors give almost identical results once saturation is properly implemented. In this way any potential ambiguity of the formalism is removed. We note also that it was  demonstrated in [15] that the explicit inclusion of the kinematical constraint again only leads to a different normalization, once saturation is implemented.
In order to demonstrate the saturation of Q s we extend in figure 20 the running coupling results to asymptotic values of Y . Once again we see the complete saturation of Q s for smaller p, but in this case one needs much larger values Y to reach the stalling region. Furthermore we see that the scaling behavior noted above in the fixed coupling case, namely that Q s (Y, p)/p saturates to a number independent than p, is no longer true when the coupling is running. Additionally, as in the previous figures, we see that the dependence on p is much weaker than in the fixed coupling case. Next, we extract Q s in the case when q is used as scale for the running coupling. As we have seen from the results for the gluon distribution, the evolution for smaller p is generally slowed down in this case, and consequently we find a slower growing saturation scale. The results are shown in figure 21. Since the evolution is generally slower, the p dependence for smaller Y is also weaker. For phenomenologically interesting Y values, however, the differences are not that large. In all the results presented here for Q s one should keep in mind that the estimates are highly optimistic, in the sense that Q s is rather large already at low Y . From phenomenological analysis performed at HERA, we do not expect to find such a large normalization of Q s at the LHC. We discuss the phenomenological relevance of our method in the next section.
Finally let us mention that, though not shown, we have also checked that the conclusions are not altered when a different set of the boundary parameters are used. In addition to the default choice ∆ = 2, c = 0.4 we have also used ∆ = 1, c = 1 and ∆ = 5, c = 0.1. The different choices lead to different normalizations of Q s but do not alter the Y or p dependence. The method is therefore rather general and robust.

Summary and outlook
In this paper we have analyzed the nonlinear effects in the CCFM equation using a boundary method. We have found that, the resulting gluon distribution exhibits similar properties as in previous BFKL/BK studies only in the region of transverse momenta lower than the scale p in CCFM. As the rapidity grows and the momenta diffuse to the higher values, the dependence on the rapidity and p becomes highly complicated, and completely novel effects are observed. As a result the diffusion of the transverse momenta to higher values is significantly suppressed due to angular ordering when the typical scales exceed p. For ultrahigh rapidities this results in the stalling of the propagation front of the solution. From the gluon distribution we have extracted the saturation scale in CCFM. We have found that due to the presence of the additional scale in CCFM, the saturation scale becomes also dependent on this scale. Also, the high rapidity behavior of the saturation scale is significantly modified in CCFM. For ultrahigh rapidities the saturation scale stops to grow completely, and this is related to the above-mentioned stalling of the gluon distribution. The smaller the value of p, the lower are the rapidities for which deviations from the BK-type behavior of the saturation scale are observed. Analytical calculations confirm the numerical analysis on the behavior of the saturation scale.
We have investigated the cases both in fixed and running coupling scenarios, and found that in the running coupling case the novel features are present, yet for the higher values of rapidities. We have also investigated the dependence of the gluon distribution on the choice of the non-Sudakov form-factors. In the case when the boundary is present there is very little sensitivity of the solutions on the different form factors, due to the fact that the nonlinear effects basically regularize the infrared region where the differences are the largest. On the other hand, in the linear case we have found that the different forms of the form factors lead to large discrepancies of the solutions.
Let us now comment on the possible phenomenological relevance and importance of the method and results presented here. The results obtained here indicate that the saturation of Q s occurs at extremely large values of rapidity. In the phenomenologically relevant case of the running coupling, the rapidities needed to observe this phenomenon are about Y = 30 − 40. Such extremely large values are not accessible at present particle colliders. One should also take into account the fact that here we only included the 1/z term of the splitting function. The other, nonsingular at z = 0 terms are important too and they are known to reduce the speed of the rapidity evolution. One also needs to properly take into account the effects of the energy momentum conservation which will further reduce the growth. Taking these effects into consideration we can conclude that it will be very difficult to observe the novel effect of saturation of Q s in the near future. One the other hand one needs to note that the exact values of the rapidities at which stalling occurs could depend on the details of the initial condition. What could be however relevant phenomenologically is the dependence of the saturation scale on the additional scale p in the regime where Q s is of order p. For example in DIS, in the usual saturation-based approaches one treats Q and Q s (x) as completely independent scales. The results of our analysis indicate that saturation scale should be correlated with Q when both scales are of the same order.
The interesting question is, however, whether or not it is urgent to introduce saturation effects for studies at LHC involving unintegrated parton densities. This question is relevant and important, as saturation affects the shape of the unintegrated distributions, even at higher k. Usually the distributions are extracted using linear evolution equations and they will be used in almost every physics analysis at the LHC. A recent study [34] has reported deviations from small-x resummed DGLAP evolution in the HERA data that might be related to saturation effects. Therefore we think that there is a need to analyze the possible saturation effects in the evolution and their impact onto the LHC observables with substantial precision.
To study in more detail the relevance of saturation, one would first have to include additional corrections to the hard emissions. Obviously the soft emissions associated with the 1/(1 − z) pole are important to include as well as the other non-singular terms in z of the splitting function. We intend to return to this issue in a future publication. The soft emissions are included in the CASCADE event generator, and as already outlined in [14,15], the boundary method can be implemented in an event generator much in the same way. A first attempt of applying a boundary method in CASCADE was reported in [35]. In this case, however, Q s was assumed a priori to be of the standard form Q s ∼ x −λ [19], and it was then used as a cut-off in k. As we have seen here, Q s in CCFM depends on rapidity and p in much more complicated way. Moreover, the power of the boundary method is that it gives a Q s consistent with the evolution.
A problem we see with the application of the boundary method in CASCADE is that the results coming from it show almost no growth at all of A from the small-x evolution. The relatively steep growth in x needed to fit data at HERA is rather provided by the x dependence of the off-shell hard matrix element and the initial condition. Obviously, if there is very little growth coming from the perturbative gluon cascade, then saturation will not be relevant.
A similar effect is observed in the resummed BFKL/DGLAP approaches [29,30,[36][37][38][39][40][41][42]. It was consistently found that the resummed gluon-gluon splitting function does indeed have the singularity at small x, but the subleading corrections cause the resummed splitting function to be very close numerically to NLO DGLAP in the region x ≥ 10 −4 . Therefore in the kinematic regime covered by HERA the evolution of the gluon distribution is controlled by the DGLAP equations. On the other hand it was found in earlier analysis [32] that when evaluating F 2 structure function from the k T factorization approach together with unintegrated gluon distribution, large enhancement in the low x region is provided by the off-shell matrix element, similarly to what was observed in CASCADE.
A related problem was also reported in [43] where some phenomenological applications of CCFM were presented. In that case, again only the hard emissions were included, but to simulate possible non-leading effects, the non-Sudakov (2.7) was modified by simply changing the integration intervals, resulting in a larger suppression coming from the form factor (the exact motivation for this change was given by checking that the resulting high energy saddle point more or less agree with the one extracted from NLL BFKL studies). The result was that the small-x growth turned out be much delayed again, and as discussed in that paper this delay might be behind the fact that CCFM predictions on forward jets seem to underestimate the growth of the jet cross section with decreasing x (for recent studies of forward jets within CASCADE see [44,45]).
On the other hand, a reason why one has to introduce a strong suppression from additional effects, which then delay the small-x growth, is because it is otherwise not possible to fit F 2 from the linear evolution. This is not so surprising since the linear evolution at small momenta gives a very strong growth which we know must be tamed by non-linear effects. Even if one only performs a fit to F 2 for moderate values of Q, because one is integrating A to smaller k, saturation effects are not completely negligible. One can therefore imagine that a procedure like that in [43], where the non-leading effects were introduced by hand, easily overestimates the suppression coming from them when saturation effects are not properly included. On the other hand the non-leading pieces are more properly included in CASCADE which has the same problems as in [43]. An event generator is, however, rather complex, including many different components, and it is not always very transparent how big the role of the first principle physics in it is. Therefore we believe it will be very valuable to include the soft emissions in a simple numerical simulation as done here, so that one can more closely study the combined effects of the different physics components in a simpler environment. This may shed light on the possible relevance of saturation for hard physics at the LHC (obviously for the possible description of the Underlying Event (UE) saturation effects such as multiple scatterings etc are very important. For a possible description of the very complex UE at the LHC one will need an entirely different approach to the problem, for some recent articles on this issue see [46,47]).