The chaotic effects in a nonlinear QCD evolution equation

The corrections of gluon fusion to the DGLAP and BFKL equations are discussed in a united partonic framework. The resulting nonlinear evolution equations are the well-known GLR-MQ-ZRS equation and a new evolution equation. Using the available saturation models as input, we find that the new evolution equation has the chaos solution with positive Lyaponov exponents in the perturbative range. We predict a new kind of shadowing caused by chaos, which blocks the QCD evolution in a critical small $x$ range. The blocking effect in the evolution equation may explain the Abelian gluon assumption and even influence our expectations to the projected Large Hadron Electron Collider (LHeC), Very Large Hadron Collider (VLHC) and the upgrade (CppC) in a circular $e^+e^-$ collider (SppC).


Introduction
The QCD evolution equation is an important part in the study of high energy physics.
As we know, the nonlinear iteration equations may have a characteristic solutionchaos, which has been observed in many natural phenomena [6]. A following question is: do the nonlinear QCD evolution equations have chaotic solution? Several years ago we have reported chaos in a new evolution equation [7], which describes the corrections of the gluon recombination to the BFKL equation at the leading logarithmic LL(1/x) approximation. The purpose of this work is to detail this discovery after a long consideration.
We begin from the proposal of the new evolution equation. Fig. 1 is a schematic program, which shows that the correlations among initial gluons modify the evolution  We will present the derivations of the above mentioned four evolution equations in a same partonic framework. For this sake, we use the Bjorken frame, where the traditional parton distributions inside a fast moving target are defined in the factorization scheme.
Note that the BFKL equation was originally derived by using the Regge langauge. In this work we take an alternative technic to re-derive the BFKL equation in Sec. 2, where the time ordered perturbation theory (TOPT) [8] is used the same as the Altarelli-Parisiderivation in the DGLAP equation [2].
The new derivation of the BFKL equation allows us conveniently to add the corrections of the gluon fusion on it according to the physical pictures in Fig. 1. We present the derivation of the new evolution equation in Sec. 3. The nonlinear part of this equation has IR divergences similar to the linear BFKL-kernel. Naturally, the similar regularization scheme as in the BFKL equation is necessary. Thus, we use the TOPT-cutting rule [4] to collect the contributions from the virtual processes in the linear and nonlinear parts of the new evolution equation. Four evolution equations at small x in Fig. 1 show their consistence. We discuss the relations among these evolution equations in Sec. 4. We find that the new evolution equation is a natural result following the DGLAP, BFKL, GLR-MQ-ZRS and BK equations.
Using the available saturation models as the input distribution, we study the numerical solutions of Eq. (3.46) in Sec. 5. The solution shows an unexpected result: the unintegrated gluon distribution function F (x, k 2 ) in Eq. (3.46) begins its smooth evolution under suppression of gluon recombination, but when x approaches a small critical value x c , F (x, k 2 ) will oscillate aperiodically in a narrow k 2 range (see Fig. 16). We find that this solution presents the chaotic characteristics. In particular, this solution of Eq.
We indicate that chaos in Eq. (3.46) origins from a serious of perturbations when k crosses over the saturation scale. The rapid oscillation in chaos in a narrow k 2 domain arises a big shadowing (Fig. 15), which blocks the QCD evolution vis three gluon vertex (Fig. 14). The chaos effects in Eq. (3.46) are discussed in Sec. 6. Chaos, which has been observed in nature, is a highlighted phenomenon in nonlinear physics. We proposed an example where chaos appears in a QCD evolution equation and it may influence the gluon distribution function, even change our expectations to the future large hadron colliders.
In this paper, sections 1-4 are the derivation of the new evolution equation; sections 5-6 present the chaos solution of this equation and its effects.
We consider the following partonic picture of the DIS process. At the lowest order, the elementary amplitude in Fig. 1a together with its conjugate amplitude constructs the DGLAP equation for gluon. However, this picture should be modified at small x due to the correlations among initial gluons. For example, a possible correction to the DGLAP-amplitudes are given in Fig. 1b, or detailed in Fig. 2. These processes imply that a scattered gluon is omitted from two correlating gluons before its radiation. We call such a correlated gluon cluster as the cold spot, which phenomenologically describes the correlation among initial partons, where the dark circle implies soft QCD-interactions.
Neglecting the irrelevant part with the evolution dynamics using the TOPT decomposition, using the TOPT-decomposition Fig. 2 can been simplified as Fig. 3, where the dashed lines are the time-ordering lines in the TOPT and "x" marks the probing place.
Note that the all lines across the time lines are on mass-shell.
The evolution kernel in QCD evolution equation is a part of a complete scattering diagram. In general, the correlated initial partons have the transverse momenta and they are off mass-shell, therefore, the k-factorization scheme is necessary. In this work we use the semi-classical Weizsäcker-Williams (W − W ) approximation [9] to realize the kfactorization scheme. The W −W approximation allows us to extract the evolution kernels and to keep all initial and final partons of the evolution kernels on their mass-shell. These diagrams lead to the real part of the BFKL equation. According to the scale-invariant parton picture of the renormalization group theory [10] the observed wave function Ψ(x 2 , k) is evolved from the initial wave functions Ψ(x 1 , p a ) and Ψ(x 1 , p b ) via the QCD interactions, i.e., where the two perturbative amplitudes corresponding to Fig. 3 are and 3) The momenta of the partons are parameterized as and The matrices of the local QCD interactions are (2.10) where the polarization vectors are and ǫ(l a ) = (0, ǫ, − ǫ · l a (x 1 − x 2 )P ), (2.13) where ǫ is the transverse polarization of the gluon in ǫ µ = (ǫ 0 , ǫ, ǫ 3 ) = (0, ǫ, 0), since the sum includes only physical transverse gluon states in the TOPT form.
Taking the LL(1/x) approximation, i.e., assuming that x 2 ≪ x 1 , one can get two similar amplitudes 14) and However, these two amplitudes really occupy different transverse configurations. This is a reason why the dipole model of the BFKL equation is derived by using the transverse coordinator-space. However, we shall show that the momentum representation still can be used to distinguish the differences between Eqs. (2.14) and (2.15).
The two parton correlation function is generally defined as (2.16) where k c and k ab are conjugate to the impact parameter and transverse scale of a cold spot.
Equation (2.16) implies the probability of finding a gluon, which carries the longitudinal momentum fraction x of a nucleon and locates inside a cold spot characterized by k c and k ab .
In this work we derive the evolution equations in the impact parameter-independent case. This approximation implies that the evolution dynamics of the partons are dominated by the internal structure of the cold spot. Thus, the evolution kernel is irrelevant to k c and we shall use (2.17) which has the following TOPT-structure Notice that all transverse momenta in Eqs. (2.4)-(2.13) are indicated relative to the mass-center of the nucleon target. However according to Eq. (2.17), the evolution variable is the relative momentum k ab , therefore, it is suitable to rewrite all momenta relative to p b in Eq. (2.2) and to p a in Eq. (2.3), respectively. Thus, we replace the transverse momenta as follows: 19) in Eq. (2.2) since 21) in Eq. (2.3). In consequence, we have and where we identify two ǫ in Eq. (2.23) since the measurements on (x 2 , k 2 a0 ) and (x 2 , k 2 0b ) are really the same event.
Equation (2.1) together with Eqs. (2.22) and (2.23) provide such a picture: a parent cold spot with the longitudinal momentum fraction x 1 and transverse momentum k ab radiates a gluon, which has the longitudinal momentum fraction x 2 and the transverse momentum k a0 (or k 0b ). It is interesting that this is a picture like the dipole model but in the full momentum space. In fact, using the Fourier transformation, one can obtain the corresponding amplitude in the dipole model [11] A BF KL (x a0 , x 0b , x 1 , ] · ǫ. (2.24) where x is the conjugate coordinator corresponding to the relative transverse momentum k.
We taking the square of the total amplitude, one can get (2.25) where the probe in the last step only picks up the contributions from Ψ(x 2 , k a0 )Ψ * (x 2 , k a0 ), we regard ∆f (x 2 , k a0 ) as the increment of the distribution f (x 1 , k ab ) when it evolves from (x 1 , k ab ) to (x 2 , k a0 ). Therefore we have According to Eq. (2.23), the evolution kernel reads as This is the real part of the BFKL equation.
The evolution kernel of the DGLAP equation has infrared (IR) singularities, which relate to the emission or absorption of quanta with zero momentum. A standard regularized method is to combine the contributions of the corresponding virtual processes.
We call a cut diagram as the virtual diagram, where one side of the cut line is a naive partonic definition without any QCD corrections. A simple calculation of the virtual diagrams was proposed via the TOPT cutting rule in [4]. Let us summarize the TOPT cutting rule as follows. When we use a probe to observe the parton distributions inside the target, we cannot control the probing position. In principle, we should sum over all cut diagrams belonging to the same time-ordered un-cut diagrams, and these diagrams have similar singular structure but may come up with opposite signs. The TOPT-cutting rule presents the simple connections among the related cut-diagrams including the real-and virtual-diagrams. The BFKL-kernel also has singularities on the transverse momentum space. Thus, we can pick up the contributions from the virtual diagrams using the TOPT-cutting rule without the complicated calculations. Using the TOPT-cutting rule, one can prove that the diagrams in Fig. 4 contribute a similar evolution kernel as the real kernel but differ by a factor −1/2 × (1/2 + 1/2).
The negative sign arises from the changes of time order in the energy denominators. The factor (1/2 + 1/2) is due to the fact that the probe "sees" only the square root of the parton distribution, which accepts the contributions of the partonic processes in a virtual diagram, and the other factor 1/2 is originated from the symmetry of the pure gluon process. Therefore, the evolution equation corresponding to Fig. 4 is (2.32) Since we calculate the contributions to ∆F (x, k a0 ), we should make the replacement b ↔ 0 in Eq. (2.32). Combining the real and virtual parts of the evolution equation, we have According to Eq. (2.18), the distribution f (x, k) in the TOPT-form contains a singular factor 1/k 4 , which arises from the off energy-shell effect in the square of the energy denominator. In order to ensure the safety using of the W − W approximation, we move this factor to the evolution kernel and use the following new definition of the unintegrated gluon distribution wherek is a unity vector on the transverse momentum space. Thus, Eq. (2.33) becomes which is consistent with a standard form of the BFKL equation.
The correlations among the initial gluons can be neglected in the dilute parton system.
In this case the contributions of the interference diagrams Figs.3c and 3d disappear. Thus, the kernel Eq. (2.30) reduces to the splitting functions in the DGLAP equation at the small x limit, i.e., Since in this case two initial gluons have the same transverse momentum, we can always take it to zero and use the collinear factorization to separate the gluon distribution. The corresponding DGLAP equation reads where the scaling restriction δ(x 2 − x B ) is included and (2.38)

The new evolution equation
We consider the evolution kernel based on Fig. 1d, which constructs a new evolution equation. Notice that the two pairs of initial gluons, which are hidden in the correlation function, for example in Fig. 5a, should be indicated as Fig. 5b.
The amplitude 3) and Fig.6: The TOPT-diagrams constructed by the elemental amplitudes in Fig. 1d. For simplicity we neglect some lines linking with l a , l b ... in 7c and 7d, since they are irrelevant to the evolution kernel. The momenta of the partons, for example, are parameterized as For example, in the t-channel and The matrices in Eqs. (3.3) and (3.4) are 14) and where d γη ⊥ = n γ n η +n η n γ −g γη , C αβγ C ρση are the triple gluon vertices and the polarization vectors are . (3.19) Thus, at small x we have 20) where one of the two factors in each term is from the approximation We use the relative transverse momenta to replace the relating momenta in Eqs. (3.5)- (3.13) and recalculate Eqs. (3.3) and (3.4). The result is where the foot-indexes of the relative transverse momenta indicate the corresponding cold spots and k(m), k(m ′ ) imply that the momenta origin from m, m ′ , respectively. Using the definitions we have We read two momenta k (p b ,pc) in Eq. (3.22) as k 0c and k b0 , respectively. On the other hand, due to momentum conservation, we have Thus, we obtain (1) The momenta p a and p d don't flow into the amplitude Eq. (3.26). Therefore, the cold spot (p a , p d ) in Fig. 6b is independent of the evolution dynamics and its distribution should be integrated as a unobservable quantity. Thus, the resulting kernel reduces to the linear BFKL kernel.
(2) The momenta p a and p d flow into the amplitude Eq. (3.26) through m and m ′ . 27) and where one can introduce In general, the momenta k a0 and k 0d in Eqs. (3.27) and (3.28) are undetermined since l a and l d in Fig. 6b are unobserved, they should be integrated out as two independent variables. Thus, the resulting evolution kernel reduce to the DGLAP-like kernel.
Obviously, the above mentioned two situations should be excluded in our resummation in order to get the leading corrections, unless we have the following restriction conditions and they imply that due to Eqs. (3.24) and (3.29). To understand Eq. (3.31), we image that before the probe interacts with the target, two overlapping cold spots have recombined into a common cold spot (p b , p c ), This is an inverse processes of the dipole splitting in the BK equation [5].
Therefore the probe always measures the recombination processes of four initial gluons originated from a same cold spot and sharing a same relative momentum.
Summing all the channels, we get the evolution kernel corresponding to Fig. 6 and the result reads In the case of decreasing gluon density, the contributions of the interference terms (Figs. 7c and 7d) disappear and Fig. 1d return to Fig. 1c. Thus, Eq. (3.32) reduces to the real part of the GLR-MQ-ZRS kernel [12] K N ew where a power suppressed factor 1/Q 2 1 has been extracted from the evolution kernel.
The correlation function G (2) is a generalization of the gluon distribution beyond the leading twist. It is usually modeled as the square of the gluon distribution. For example, (3.35) where R N is the correlation scale of the gluons in the nucleon. The definition (3.35) is a phenomenological model, which contains an arbitral normalization constant. However, this constant will be determined through the value of R N by using the experimental data.
Where we meet very complicated calculations about the interference-and corresponding virtual amplitudes. However, the TOPT-cutting rule shows that the above mentioned amplitudes correspond to a similar recombination kernel except the numerical factor and the different kinematic regions [4].
Another key problem is that we meet various multi-gluon correlation functions, in which the cut line cuts off the nonperturbative matrix with different ways. Fortunately, Jaffe has shown that these correlation functions on the light-cone has the same form in the DIS processes [13]. The Jaffe-cutting rule was broadly used in the study of the high twist processes. The TOPT provides a straightforward explanation about the Jaffe-cutting rule: since all backward propagators are absorbed into the nonperturbative correlation functions, the partons correlating two initial gluons inside the nonperturbative matrix are on mass-shell. Therefore, the correlation functions with cuts at different places are the same. Thus, the Jaffe-cutting rule can be included in our TOPT-cutting rule. Combining the DGLAP dynamics at small x, the GLR-MQ-ZRS equation reads where we take the same parameter R N as in Eq. (3.35) since the relation (2.28) is irrelevant to R N . Using the evolution kernel (3.32), we write Now let us discuss the contributions from the virtual diagrams. According to the standard regularization schema, the TOPT-cutting rule shows that the diagrams in Fig.   7 have a similar evolution kernel as that in Fig. 6 but with the different kinematical variables and differ from a simple numerical factor.
The processes in Figs. 6 and 7 contributes the net positive antishadowing effect. The negative shadowing effect is really originated from the interference processes, two of them are shown in Fig. 8. Here the contributions from the corresponding virtual processes are also necessary (see Fig. 9). The TOPT-cutting rule shows that the processes in Figs. 8 and 9 also have a similar evolution kernel.  Up to now we have separately established the relations of the evolution kernels between the real and virtual diagrams in the 4-partons-to-4 partons (4 → 4) amplitude and the 3-partons-to-5-partons (3 → 5) amplitude, respectively. In the next step we will show that the relationship between the above mentioned two kinds of virtual diagrams will link up all the four evolution kernels. According to Eq. (3.20), the resulting amplitudes are irrelevant to the transverse momenta of the initial gluons at x 2 ≪ x 1 . Thus, we use the relations shown in Fig. 10, which are derived in the collinear factorization schema [4] to reveal that the two kinds of virtual diagrams differ only from a minus sign, which is from an energy deficit between the two dashed lines in Fig. 9: because both the momenta k b0 and k 0c are indicated by k in the mass-center of the nucleon target, we have on the left-hand side of Fig. 10, where x m < x l , (x m and x l are the longitudinal momentum fractions in the momenta m and l, respectively); and on the right-hand side of Fig. 10, where x m > x l . In consequence, we finally link up all evolution kernels and obtain the following equa- where we assume that the Jaffe-cutting rule is still holden in the k-factorization scheme Thus, the correlation function can be cut and we can use the same correlation function in the real, virtual, and interference processes. From Eq. (3.41) we have Similar to Eq. (2.34) we note that Combining it with the linear BFKL equation, we finally obtain a complete evolution .

(3.46)
Comparing with the GLR-MQ-ZRS equation (3.36), the contributions of the virtual diagrams can't be canceled in Eq. (3.46) and they are necessary for IR safety.

Unity of the QCD evolution equations
It is a surprise that Eq. (3.46) can be "directly" written by using an analogy with the DGLAP, BFKL and GLR-MQ-ZRS equations. For this sake, we summarize the four evolution equations at small x as follows. The DGLAP equation ( (see Fig. 12b); The GLR-MQ-ZRS equation (3.36) (see Fig. 12c); The equation (3.46) Fig. 12d).
One can find the following interesting relations among these equations: The DGLAP and BFKL equations have the same evolution dynamics (i.e., the gluon splitting), where we have the following analogy between the real parts of Eqs. (4.1) and (4.2): The nonlinear parts of the GLR-MQ-ZRS and Eq. (3.46) also have the same evolution dynamics (i.e., the gluon recombination), they have similar relationships like Eqs. (4.5) and (4.6): G(x, k 2 ) ↔ F (x, k bc ), (4.8) and an extra relation for the power suppression factor 1 k 2 ↔ The BK equation [5] is generally considered as a typical nonlinear correction to the BFKL equation at the LL(1/x) approximation. We discuss the relation of Eq. (3.46) with the BK equation. The BK equation is usually written by using the scattering amplitude (4.14) The nonlinear evolution kernel in the BK equation is regularized by the connecting amplitude N(x bc , x)N(x c0 , x) rather than using the virtual diagrams. Using (4.15) and the definition (4.16) one can obtain the BK equation in the momentum space Since the measured unintegrated gluon distribution F (x, k 2 ) is irrelevant to the azimuthal angle φ (see Eq. 2.38), after azimuthal integration we have where we use Eq. (3.33), i.e., (4.20) The nonlinear evolution kernel in Eq. (4.19) is essentially the GLR-MQ-ZRS-kernel [12] and it collects only the k 2 -ordered corrections. That is, k 2 bc are ordered in [k 2 min , k 2 b0 ].
As an approximation, we only keep the last step evolution, i.e., we set k 2 bc = k 2 b0 and call it as the one step evolution approximation. Insert the dimensionless function δ(1 − k 2 b0 /k 2 bc ) into Eq. (4.19), (4.21) which leads to the BK-like equation (4.17) where we take We will focus on the behavior of the solutions near the saturation range, where we estimate that F (x/2, k 2 ) ≃ F (x, k 2 ) in Eq. (3.46) due to the strong shadowing effect.
Thus, Eq. (3.46) reduces to The solutions of Eq.(5.2) depend on the strength of the nonlinear terms, which include the model-dependent assumptions in Eqs. (3.35), (3.37) and a free parameter R N . To reduce the uncertainty, the value of R N = 4GeV −1 with the assumption (3.37) is independently fixed by fitting the available experimental data about the proton structure function using the GLR-MQ-ZRS equation in Ref. [15].
The solutions of Eq. (5.2) need the knowledge of the gluon distribution with all k 2 at a starting x 0 . A major difficulty is the treatment of the infrared region, k 2 < k 2 0 (k 2 0 ∼ 1 GeV 2 ). The BFKL evolution leads to diffusion of the starting k-distribution both to larger and to smaller values of k. However, the perturbative BFKL-growth of F (x, k 2 ) toward smaller k 2 is not expected to be valid when the gluon momenta enter the nonperturbative region. The common feature of nonperturbative modifications of the infrared region is that the solution F (x, k 2 ) vanishes as k 2 → 0. The reasons, for example, are the requirement of gauge invariance [16], the colour neutrality of the probed proton [17], and the absence of the valence gluons in a static proton [15]. Therefore, the increasing distribution F (x, k 2 ) should be saturated at k 2 < Q 2 s (x), Q s (x) is called as the saturation scale. For example, with the color-dipole approach Golec-Biernat and Wusthoff (GBW) [18] used the inclusive and diffractive scattering data and obtained where σ 0 = 29.12mb, x 0 = 0.4 × 10 −4 , λ = 0.277, R 0 (x) = (x/x 0 ) λ/2 /Q s and Q s = 1GeV .
Note F ≡ F/k 2 in Eq. (2.38). The parameter α s is fixed as α s = 0.2. The GBW model gives a description of F near the saturation scale, although it lacks the QCD evolution.
We draw F input (x 0 , k 2 ) = k 2 F GBW (x 0 , k 2 ) in Fig. 13 (the solid curve). In the calculations we divide the evolution region into two parts: region(A) 0 to Q 2 s and region(B) Q 2 s to ∞. In region(B) the QCD evolution equation is taken to evolute and in region(A) the nonperturbative part of F (x, k 2 ) is identified as (5.4) where the parameter C keeps the connection between two parts.
The Runge-Kutta method is used to compute Eq. (5.2). Note that F (x,   The sudden change of a solution is an interesting phenomenon in nonlinear evolution system, in particular, this behavior perhaps relates to chaos. An important character of chaos is that the solution is sensitively relevant to the initial conditions. For this sake, we study the solutions of different input conditions. We compute a similar solution as Fig.   16 but the starting point is moved a little from x 0 = 0.4 × 10 −4 to x 0 = 0.35 × 10 −4 . The results in Fig. 17 show that the oscillation structure of F (x, k 2 ) ∼ k 2 is sensitive to the starting point of the evolution, although the global behaviors of the curves are similar.
We change the input distribution to 1.01×Eq. (5.4) and compare these results in Fig.   18. One can find the obvious difference in the oscillation structure. The improvement of the precision in the computation may aggravate the chaotic oscillations since the increasing samples perturb the distributions at every step in the evolution.
In contrast, if the above mentioned oscillations are arisen from the calculation errors, such oscillations will disappear with the increasing precision. In Fig. 19 we present the curve with a same input as in Fig. 16 but with double calculating precision. One can find that the oscillations are aggravated with increasing precision.
The above aperiodic oscillation is sensitive to the initial conditions. Especially, the oscillation will be enhanced with the increase of the numerical calculation precision. These features are universally observed in many chaos phenomena.
A standard criterion of chaos is that the system has the positive Lyapunov exponents, which indicates a strong sensitivity to small changes in the initial conditions [6]. We regard y = ln 1/x as 'time' and calculate the Lyapunov exponents λ(k 2 ) in a finite region, where the distribution oscillation is obvious. We divide equally the above mentioned y-region into n parts with y 1 , y 2 .., y n+1 and τ = (y n+1 − y 1 )/n. Assuming that the distribution evolves to y 1 from y 0 = ln 1/x 0 and results F (y 1 , k). Corresponding to a given value F (y 1 , k) at (y 1 , k), we perturb it to F (y 1 , k) + ∆ with ∆ ≪ 1. Then we continue the evolutions from F (y 1 , k) and F (y 1 , k) + ∆ to y 2 from y 1 respectively, and denote the resulting distributions as F (y 2 , k) andF (y 2 , k). Making the difference ∆ 2 = |F (y 2 , k) − F (y 2 , k)|. In the following step, we repeat the perturbation F (y 2 , k) → F (y 2 , k) + ∆ and let the next evolutions from F (y 2 , k) and F (y 2 , k) + ∆ from y 2 to y 3 respectively and get the results ∆ 3 = |F (y 3 , k) − F (y 3 , k)|...... (see Fig. 20). The Lyapunov exponents for the image from y to F (y, k) are defined as The Lyapunov exponents of the gluon distribution in Eq. (5.2) with the input Eq. (5.4) are presented in Fig. 21. For comparison, we give the Lyapunov exponents of the BFKL and BK-like equations. The positive values of the Lyapunov exponents clearly show that the oscillation of F (x, k)∼k 2 is chaos of Eq. (5.2). Therefore, we conclude that chaos in Eq. (5.2) blocks the QCD evolution of the gluon distribution. Fig.20: Schematic programs to calculate the Lyapunov exponents of the evolution equations.

Discussions
The exact value of x c depends on the initial conditions, which have some uncertainties, however, the fact of chaos is irrelevant to the detailed dynamics, provided an essential change of the k 2 -dependence of F (x, k 2 ) when the evolution transfers from perturbative to nonperturbative ranges. For example, an alternative saturation model [19] assumes with k 2 a = 1GeV 2 to replace Eq. (5.4) (for x = x 0 , see the dashed curve in Fig. 13).
The chaos solutions still exist in Fig. 22 where x c ∼ 1.7 × 10 −7 . The reason of the chaos solution in Eq. (5.2) is that this equation contains the following regularized kernels in the linear terms and  We discuss qualitatively the azimuthal angle (φ)-dependent case. Equation (5.2) be- ..... , (6.5) where (......) are the non-singular parts. One can find that the equation contains the similar regularized forms like Eqs. (6.3) and (6.4): (6.6) in the linear terms and in the nonlinear terms. As we have emphasize that both the chaos and the blocking effect origin from such kind of regularized kernels. Therefore, we consider that our results are still hold for the azimuthal angel-dependent solutions.
The equation (3.46) is based on the leading QCD corrections, where the higher order corrections are neglected. An important questions is: will the chaos effects in the new evolution equation disappear after considering higher order corrections? We have known that chaos in the MD-BFKL equation origins from the singularity of the nonlinear evolution kernel. From the experiences in the study of the BFKL equation, higher order QCD corrections can not remove the singularities at the lower order approximation [20].
In particular, the virtual cut diagrams always exit in any higher order corrections to the BFKL equation. The regularization similar to Eq. (6.4) is necessary. Besides, the chaotic behavior cannot be destroyed by arbitrarily small perturbations of the system parameters.
Therefore, we expect that chaos still exists in Eq. (3.46) even considering the higher order corrections.
The solution F (x, k 2 ) of Eq. (5.2) becomes zero can not be simply explained as the gluon disappearance at x ≤ x c . Although the three gluons vertex stops working at x < x c , the gluons still can evolve similar to the Abliean photons in a thin parton system.
In a quark confinement mechanism, the dual-superconductor picture was suggested by Refs. [21][22][23], where an assumption of Abelian dominance seems to be significant to confinement. The Abelian dominance means that only the diagonal gluon component in the confinement mechanism. The distributions of the non-Abelian gluons collapse at x < x c , the contributions of the Abelian gluons appear. One can image that the Abelian gluons dominate the soft gluons. Thus, chaos in Eq. (3.46) provides a dynamical mechanism for separating the Abelian gluons.
The blocking effect in the QCD evolution will suppress the new particle events in an ultra high energy hadron collision. Although we have not exactly predicted the energy scale x c which corresponding to the blocking effect, the chaos solutions in Eq. (3.46) should arise our attention when considering the future large hadron collider. In particular, the nonlinear coefficients in the evolution equation will be enhanced by a factor [1+0.21(A 1/3 − 1)] in the nuclear target since the correlations of gluons among different bound nucleons [24], this will increase the value of x c into the observable range of the projected Large Hadron Electron Collider (LHeC) [25], Very Large Hadron Collider (100TeV VLHC) [26] and the upgrade (CepC, CppC) in a circular e + e − collider (SppC) [27]. Figure 23 presents the nuclear A-dependence of x c . We will detail them elsewhere. Although the position of chaos is undetermined due to the value of x c sensitively dependent on the input conditions, the existence of chaos in the QCD evolution equation may change our expectation to the future large hadron collider plans.