Nonlinear GLR-MQ evolution equation and Q^2-evolution of gluon distribution function

In this paper we have solved the nonlinear Gribov-Levin-Ryskin-Mueller-Qiu (GLR-MQ) evolution equation for gluon distribution function G(x,Q^2) and studied the effects of the nonlinear GLR-MQ corrections to the Leading Order (LO) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations. Here we incorporate a Regge like behaviour of gluon distribution function to obtain the solution of GLR-MQ evolution equation. We have also investigated the Q^2-dependence of gluon distribution function from the solution of GLR-MQ evolution equation. Moreover it is interesting to observe from our results that nonlinearities increase with decreasing correlation radius (R) between two interacting gluons. Results also confirm that the steep behavior of gluon distribution function is observed at R=5 GeV^{-1}, whereas it is lowered at R=2 GeV^{-1} with decreasing x as Q^2 increases. In this work we have also checked the sensitivity of \lambda_G in our calculations. Our computed results are compared with those obtained by the global DGLAP fits to the parton distribution functions viz. GRV, MRST, MSTW and with the EHKQS model.


Introduction
The small-x, where x is the Bjorken scaling variable, behavior of quark and gluon densities is one of the challenging problems of quantum chromodynamics (QCD). The most important phenomena in the region of small-x which determine the physical picture of the parton (quark and gluon) evolution or cascade, are the increase of the parton density at x → 0, the growth of the mean transverse momentum of a parton inside the parton cascade at smallx, and the saturation of the parton density [1]. The parton distributions in hadrons play a key role in understanding the standard model processes and in the predictions for such processes at accelerators. Therefore the determination of parton densities or more importantly the gluon density in the small-x region is particularly interesting because here gluons are expected to dominate the proton structure function. The study of gluon distribution function is also very important because it is the basic ingredient in the calculations of different high-energy hadronic processes like mini jet production, growth of total hadronic processes etc. Moreover precise knowledge of the gluon distribution at small-x is essential for reliable predictions of important p-p, p-A and A-A processes studied at the relativistic heavy-ion collider (RHIC) [2] and at CERN ′ s large hadron collider (LHC) [3]. Knowledge of gluon density is also important for the computation of inclusive cross-sections of hard, collinearly factorizable, processes in hadronic collisions.
The most precise determinations of the gluon momentum distribution in the proton can be obtained from a measurement of the deep inelastic scattering (DIS) proton structure function F 2 (x, Q 2 ) and its scaling violation.
The measurement of the proton structure function F 2 (x, Q 2 ) by H1 [4] and ZEUS [5] at HERA over a broad kinematic region has made it possible to know about the gluon in the formerly unexplored region of x and Q 2 where, Q 2 is the virtuality of the exchanged virtual photon. This method is however indirect because F 2 (x, Q 2 ) at low values of x actually probes the sea quark distributions which are related via the QCD evolution equations to the gluon distribution. More direct determinations of the gluon distribution can be obtained by reconstruction of the kinematics of the interacting partons from the measurement of the hadronic final state in gluon induced processes. They are subject to different systematic effects and provide an substantive test of perturbative QCD. Direct gluon density determinations have been carried out using events with J/ψ mesons in the final state [6] and dijet events [7].
In perturbative QCD, the high-Q 2 behavior of DIS is given by the linear Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [8]. The number density of gluons, G(x, Q 2 ), and quarks, q(x, Q 2 ), in a hadron can be evaluated at large-Q 2 by solving the linear DGLAP equation to calculate the emission of additional quarks and gluons compared to some given initial distributions. The results are adjusted to fit the experimental data (mainly at small-x) for the proton structure function F 2 (x, Q 2 ) measured in DIS, over a large domain of values of x and Q 2 by adjusting the parameters in the initial parton distributions. Consecuently, the approximate analytical solutions of DGLAP evolution equations have been reported in recent years with significant phenomenological success [9][10][11].
DGLAP equation predicts a sharp growth of the gluon distribution function as x grows smaller which is also clearly observed in DIS experiments at HERA. This sharp growth of the gluon distribution function will have to eventually slow down in order to not violate unitarity bound [12] on physical cross sections. It is a known fact that the hadronic cross sections comply with the Froissart bound [12] which derives from the general assumptions of the analyticity and unitarity of the scattering amplitude. The Froissart bound indicates that the total cross section does not grow faster than the logarithm squared of the energy i.e., σ total = π m 2 π (lns) 2 , where, m π is the scale of the range of the strong force [13]. Gluon recombination is commonly believed to provide the mechanism responsible for the unitarization of the cross section at high energies or a possible saturation of the gluon distribution function at small-x. In other words, the number of gluons at small-x will be so large that they will spatially overlap and therefore, gluon recombination will be as important as gluon splitting. In the derivation of the linear DGLAP The proton structure function F 2 (x, Q 2 ) has been measured down to x ∼ 10 −5 but still in the perturbatively accessible region by the H1 Collaboration at HERA [4]. These data have been included in the recent global analyzes by the MRST [14] and CTEQ [15] collaborations. DGLAP evolution equations can describe the available experimental data quite well in a fairly broad range of x and Q 2 with appropriate parameterizations. But DGLAP approach cannot provide a good description while trying to fit the H1 data simultaneously in the region of large-Q 2 (Q 2 > 4 GeV 2 ) and in the region of small-Q 2 (1.5 GeV 2 < Q 2 < 4 GeV 2 ) [4,16]. This implies that towards smaller values of x and (or) Q 2 (but still Q 2 ≥ Λ 2 , Λ being the QCD cut off papameter) it is possible to observe gluon recombination effects which lead to nonlinear power corrections to the DGLAP equations. These nonlinear terms lower the growth of the gluon distribution in this kinematic region where α S is still small but the density of partons becomes very large.
Therefore, the corrections of the higher order QCD effects, which suppress or shadow the growth of parton densities, become a center of intensive study in the last few years.
Gribov, Levin, Ryskin, Mueller and Qiu (GLRMQ) performed a detailed study of this region in their pioneering papers and they suggested that these shadowing corrections could be expressed in a new evolution equation known as the GLR-MQ equation [17][18]. This equation, involves a new quantity, G 2 (x, Q 2 ), the two-gluon distribution per unit area of the hadron. The main features of this equation are that it predicts a saturation of the gluon distribution at very small-x, it predicts a critical line separating the perturbative regime from the saturation regime and it is only valid in the border of this critical line [16,19]. It is an amazing property of GLR-MQ equation is that it introduces a characteristic momentum scale Q 2 s , which is a measure of the density of the saturated gluons. It grows rapidly with energy, and it is proportional to 1/x λ with λ = 0.2 [20]. Gribov, Levin and Ryskin first suggest a nonlinear evolution equation, in which the evolution kernels, which they called as the gluon recombination functions, are constructed by the fan diagrams [17]. Later Mueller and Qiu calculated the gluon recombination functions at the double leading logarithmic approximation (DLLA) in a covariant perturbation framework [18].
The GLR-MQ equation is broadly regarded as a key link from perturbation region to non-perturbation region. There has been much work inspired by the approach of GLR-MQ which show that gluon recombination leads to saturation of gluon density at small-x [21][22]. The predictions of the GLR-MQ equation for the gluon saturation scale were studied in Ref. [16]. A new evolution equation named as modified DGLAP equation is derived by Zhu and Ruan [23] where the applications of the AGK (Abramovsky-Gribov-Kancheli) cutting rule [24] in the GLR-MQ equation was argued in a more general consideration. Here the Feynman diagrams are summed in a quantum field theory framework instead of the AGK cutting rule. In Ref. [25] parton distribution functions in the small-x region are numerically predicted by using a modified DGLAP equation with the GRV-like input distributions.
Moreover, some studies of the GLR-MQ terms in the framework of extracting the PDFs of the free proton can be found in Ref. [26]. Also other nonlinear evolution equations relevant at high gluon densities have been derived in the recent years, and the structure functions from DIS have been analyzed in the context of saturation models [27][28][29].
The solution of the GLR-MQ equation is particularly important for understanding the nonlinear effects of gluon-gluon fusion due to the high gluon density at small enough x. The solution of nonlinear evolution equations also provides the determination of the saturation momentum that incorporates physics in addition to that of the linear evolution equations commonly used to fit DIS data. Various studies on the solutions and viable generalizations of the GLR-MQ equation have been done in great detail in the last few years [30][31][32]. In the present work we intend to obtain a solution of the nonlinear GLR-MQ evolution equation for the calculation of gluon distribution function in leading order. This paper addresses interesting questions about validity of the well known Regee like parametrization in the region of moderate virtuality of photon. Here we have also calculated the Q 2 -evolution of gluon distribution function and the results are compared with the predictions of different paramerizations like GRV1998LO [33], MRST2001LO [14], MSTW2008LO [34] and EHKQS molel [16]. Finally, we present our conclusions.

Theory
The GLR-MQ equation is based on two processes in the parton cascade: the emission induced by the QCD vertex g→g + g with a probability which is proportional to α S ρ and the annihilation of a gluon by the same vertex g + g→g with a probability which is proportional to is the density of the gluon in the transverse plane, πR 2 is the target area, and R is the correlation radius between two interacting gluons. Normally, this radius should be smaller than the radius of a hadron. It is worthwhile to mention that R is non-perturbative in nature and therefore all physics that happens at distance scales larger than R is non-perturbative [30]. Here, r is the size of the parton (gluon) produced in the annihilation process. For DIS r∝ 1 Q 2 . Clearly, at x ∼ 1 only the production of new partons (emission) is essential because ρ≪1 , but at x→0 the value of ρ becomes so large that the annihilation of partons becomes important.
To take interaction and recombination of partons (mainly gluons) into account, a small parameter is introduced which enables us to estimate the accuracy of the calculation, given as, In terms of gluon distribution function this equation can be expressed as which is named as the GLR-MQ evolution equation. The factor γ is found to be γ = 81 16 for N c = 3, as calculated by Mueller and Qiu [18]. Now to study the Q 2 -evolution of gluon distribution function, we can rewrite Eq. (3) in a convenient form [31] ∂G(x, Q 2 ) where the first term in the r.h.s. is the usual linear DGLAP term in the double leading logarithmic approximation and the second term is nonlinear in gluon density.
Here, the representation for the gluon distribution G(x, To simplify our calculations we consider a variable t such that t = ln Q 2 Λ 2 , where Λ is the QCD cut off parameter. Then Eq. (4) becomes As gluons are the dominant parton at small-x, therefore, ignoring the quark contribution to the gluon distribution function the first term in the r.h.s. of Eq. (5) can be expressed as [35] ∂G(x, t) ∂t DGLAP = α S (t) 11 12 The strong coupling constant α S (t) in leading order has the form [35] α S (t) = 4π where, is the one-loop corrections to the QCD β-function and N f being the number of quark flavor. Here we consider N c = 3, and T f = 1 2 N f and N f = 4.
At small-x, the behavior of structure functions is well explained in terms of Regge-like behavior [36,37]. The small-x behaviour of structure functions for fixed Q 2 reflects the high-energy behavior of the total cross section with increasing total CM energy squared s 2 , since s 2 = Q 2 ( 1 x − 1) [38]. The Regge pole exchange picture [37] would therefore appear quite appropriate for the theoretical description of this behaviour. The Regge behavior of the sea-quark and antiquark distribution for small-x is given by q sea (x) ∼ x −α P corresponding to a pomeron exchange with an intercept of α P = 1. But the valence-quark distribution for small x given by q val (x) ∼ x −α R corresponds to a reggeon exchange with an intercept of α R = 0.5. The x dependence of the parton densities is often assumed at moderate Q 2 and thus the leading order calculations in ln(1/x) with fixed α S predict a steep power-law behavior of xg(x, Q 2 ) ∼ x −λ G , where λ G = (3α S /π)4ln2 ≃ 0.5 for α S ≃ 0.2, as appropriate for Q 2 ∼ 4GeV 2 .
Moreover the Regge theory provides extremely naive and frugal parameterization of all total cross sections [39,40]. It is suggested in Refs. [41,42] that is feasible to use Regge theory for the study of DGLAP evolution equations. The tactics for the determination of the gluon distribution function with the nonlinear correction is also based on the Regge-like behavior [43].
The Regee behavior is believed to be valid at small-x and at some intermediate Q 2 , where Q 2 must be small, but not so small that α S (Q 2 ) is too large [44,45]. Moreover, as discussed in [40] the Regge theory is supposed to be applicable if W 2 is much greater than all the other variables and so, models based upon this idea have been successful in describing the DIS cross-section when x is small enough (x < 0.01), whatever be the value of Q 2 . [20,46]. Therefore, to solve the GLR-MQ equation, we consider a simple form of Regee like behavior for the determination of the gluon distribution function at small-x given as and where M(t) is a function of t and λ G is the Regge intercept for gluon distribution function. This form of Regge behaviour is well supported by the work of the authors in Refs. [40,47,48]. According to Regge theory, the high energy i.e. small-x behaviour of both gluons and sea quarks are controlled by the same singularity factor in the complex angular momentum plane [37]. Moreover, as the values of Regge intercepts for all the spin-independent singlet, non-singlet and gluon structure functions should be close to 0.5 in quite a broad range of small-x [48], we would also expect that our theoretical results are best fitted to those of the experimental data and parameterization at λ G ≈ 0.5, where λ G is the Regge intercepts for gluon distribution function.
Substituting Eqs. (6), (10) and (11) in Eq. (5) we get Performing the integrations and rearranging the terms, Eq. (12) takes the with, and Eq. (13) is a partial differential equation which can be solved as where Γ is the incomplete gamma function and C is a constant. Although Regge behavior is not in agreement with the double-leading-logarithmic solution, namely, G(x, t) ∝ exp[Cln(t)ln(1/x)] 1/2 , but, the range where x is small and Q 2 is not very large is actually the Regge regime. Accordingly However, in the region where t is not very large the corrections for the nonlinear term in Eq. (16) can not be neglected and therefore Eq. (16) does not reduce to Eq. (17). Thus we can expect that the solution given by Eq. (16) is only valid in the region of small-x and intermediate values of Q 2 (or t).
Now, to determine the Q 2 -dependence of G(x, Q 2 ), we apply initial con- from which we obtain the value of the constant C as, From this equation the constant C can be evaluated by considering an appropriate input distribution G(x, t 0 ) at a given value of Q 2 0 . Now substituting C from Eq. (19) in Eq. (16) we obtain the Q 2 -evolution of gluon distribution function for fixed x in leading order as .
Thus we have obtained an expression for the Q 2 -evolution of gluon distribution function G(x, t) in leading order by solving the nonlinear GLR-MQ evolution equation semi-numerically. From the final expression given by Eq.
(20) we can easily calculate the Q 2 -evolution of G(x, Q 2 ) for a particular value of x by taking an appropiate input distribution at a given value of Q 2 0 .

Result and discussion
In this paper we have solved the nonlinear GLR-MQ evolution equation in order to determine the Q 2 -dependence of gluon distribution function include the new precise data on DIS from HERA together with constraints from hard scattering data. MSTW2008 presented an updated parton distribution functions determined from global analysis of hard-scattering data within the standard framework of leading-twist fixed-order collinear factorisation in the MS scheme. These parton distributions supersede the previously available MRST sets and can be used for the first LHC data taking and for the associated theoretical calculations. In EHKQS model [16] the effects of the first nonlinear corrections to the DGLAP evolution equations are studied by using the recent HERA data for the structure function F 2 (x, Q 2 ) of the free proton and the parton distributions from CTEQ5L [14] and CTEQ6L [14] as a baseline. By requiring a good fit to the H1 data, they determine initial parton distributions at Q 2 0 = 1.4 GeV 2 for the nonlinear scale evolution. In Ref. [16] it is shown that the nonlinear corrections enhance the agreement with the F 2 (x, Q 2 ) data in the region of x ∼ 3x10 −5 and Q 2 ∼ 1.5 GeV 2 . for x = 10 −2 , 10 −3 , 10 −4 and 10 −5 respectively. In all graphs the input distribution G(x, t 0 ) at a given value of Q 2 0 is taken from the GRV1998LO to test the Q 2 -evolution of G(x, Q 2 ). In our analysis we consider the kinematic range 2 GeV 2 ≤ Q 2 ≤ 20 GeV 2 , where we expect our solution to be valid.
The average value of Λ in our calcultaion is taken to be 0.192 GeV. It is observed from the figures that our results show almost similar behaviour with those obtained from different global parametrizations and also with EHKQS model.
We have also investigated the effect of nonlinearity in our results for R = 2 GeV −1 and R = 5 GeV −1 respectively. For this analysis our computed values of G(x, Q 2 ) for R = 2 GeV −1 and R = 5 GeV −1 respectively from Eq. (20) are plotted against Q 2 in Fig. 2(a) for x = 10 −2 , 10 −3 , 10 −4 and 10 −5 respectively. Here, the input distribution is taken from MSTW2008 global parametrization for a given value of Q 2 0 . We have also performed an analysis to check the sensitivity of the free parameter λ G in our results. Fig. 2(b) represents the results for the Q 2 -dependence of G(x, Q 2 ) obtained from the solution of nonlinear GLR-MQ equation given by Eq. (20) for three different values of λ G and we observed that results are very sensitive to λ G as x decreases.

Conclusion
We solve the nonlinear GLR-MQ evolution equation by considering the Refs. [49 -51]. We are also interested to obtain a solution of the nonlinear GLR-MQ equation in this form and planning to produce in a future paper.
We observe that the gluon distribution function increases with increasing Q 2 as usual which is in agreement with perturbative QCD fits at small-x, but with the inclusion of the nonlinear terms, Q 2 -evolution of G(x, Q 2 ) is slowed down relative to DGLAP gluon distribution. For the gluon distribution the nonlinear effects are found to play an increasingly important role at x ≤ 10 −3 .
The nonlinearities, however, vanish rapidly at larger values of x. It is also interesting to observe that nonlinearity increases with decreasing value of R as expected. The differences between the data at R = 2 GeV −1 and at R = 5 GeV −1 increase as x decreases which is very clear from Fig. 2(a).
Results also confirm that the steep behavior of gluon distribution function is observed at R = 5 GeV −1 , whereas it is lowered at R = 2 GeV −1 with decreasing x as Q 2 increases. We have also investigated the sensitivity of λ G in our calculations and found that results are highly sensitive to λ G as x goes on decreasing.