Nonlinear corrections for the nuclear gluon distribution in eA processes

An analytical study with respect to the nonlinear corrections for the nuclear gluon distribution function in the next-to-leading order approximation at small x is presented. We consider the nonlinear corrections to the nuclear gluon distribution functions at low values of x and using the parametrization and the nuclear modification factors obtained from the Khanpour-Soleymaninia-Atashbar-Spiesberger-Guzey model. The CT18 gluon distribution is used for the baseline proton gluon density at . We discuss the behavior of the gluon densities in the next-to-leading order and the next-to-next-to-leading order approximations at the initial scale , as well as the modifications due to the nonlinear corrections. We find that the QCD nonlinear corrections are more significant for the next-to-leading order accuracy than the next-to-next-to-leading order for light and heavy nuclei. The results of the nonlinear GLR-MQ evolution equation are similar to those obtained with the Rausch-Guzey-Klasen gluon upward and downward evolutions within the uncertainties. The magnitude of the gluon distribution with the nonlinear corrections increases with a decrease in x and an increase in atomic number A.


Introduction
The dynamics of parton interactions and the partonic structure of nuclei are prime subjects of research for both particle and nuclear physics.The formation of quark-gluon plasma inside nuclei is explored during the very first fractions of fm/c in high-energy nuclear collisions.This probe is due to a large momentum (or mass) Q≫Λ QCD scale, which is the main motivation for studying nuclear parton distributions.According to the knowledge of the parton distribution functions (PDFs) of free nucleons which comes from the measurements of deeply inelastic scattering (DIS) in lepton-nucleon (lN ) collisions, the program of extracting nuclear PDFs (nPDFs) also relies on the DIS data [1][2][3].The HERA data for the free proton reached x∼10 −5 in perturbative values of Q 2 , while the DIS-measurements for nuclear targets are bound to severely higher momentum fractions, x 10 −2 .In Ref. [4], the authors studied the prospects for constraining the nuclear parton distribution functions by small-x deep inelastic scattering at the Large Hadron Electron Collider (LHeC) [5] where its extension of the kinematic covers 4 orders of magnitude in DIS.The effect of high-precision DIS-measurements at the LHeC in Ref. [4] is illustrated by the ratio of the reduced, inclusive DIS cross-sections, σ A reduced (x, Q 2 )/σ p reduced (x, Q 2 ), where where x, y and Q 2 are the standard DIS variable, and A is the number of nucleons in a nuclear target.The LHeC promises the equivalent of 1 fb −1 of luminosity for ePb collisions at LH(e)C energies.With its large Q 2 and 1/x range nuclear shadowing can be measured very precisely.At high energies, nuclear shadowing is controlled by coherence effects.Namely, shadowing is possible only if the coherence time exceeds the mean inter-nucleon spacing in nuclei and shadowing saturates if the coherence time substantially exceeds the nuclear radius [6][7][8].Nuclear shadowing at small x (i.e., x 0.1) is experimentally well studied by NMC [9].Experiments at CERN and Fermilab focus especially on the region of small values of the Bjorken variable x and show a systematic reduction of the nuclear structure function F A 2 (x, Q 2 )/A with respect to the free proton structure function F p 2 (x, Q 2 ).This phenomenon is known as nuclear shadowing effect and is associated to the modification of the target parton distributions so that xf A i (x, Q 2 ) < Axf p i (x, Q 2 ), f i = q, g, .. [10].The relation of the bound-proton PDFs with respect to free-proton PDFs f p i is often expressed in terms of the nuclear modification factors . For a nucleus A with Z protons and N = A − Z neutrons, an average PDF is obtained as where f p/A i are the PDFs of a bound proton and the neutron contents f n/A i are obtained from f p/A i via isospin symmetry [11][12][13][14].As revealed by DIS experiments, the bound nucleon PDFs are not the same as those of a free proton, but are modified in a nontrivial way and obey the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution [15][16][17][18], which describes how the PDFs depend on the factorization scale with splitting functions P ij governing the scale evolution.For evolution of the PDFs due to the evolution equation (i.e., Eq.( 3)), a non-perturbative input at some initial scale is required to obtain a PDF set.The baseline parton distributions of a proton are parametrized in the following formal form where the coefficients α 1 and α 2 control the asymptotic behavior of f i (x, Q 2 0 ) in the limits x→0 and 1, and P i is a sum of Bernstein polynomials dependent on y = f (x) which is very flexible across the whole interval 0 < x < 1.For PDFs of nuclei, an additional dependence on the atomic mass A is required [12,14,19].In Ref. [20], the authors are discussed the nuclear cross section in terms of nuclear volume and surface contributions Therefore, the cross section per nucleon is assumed to be proportional to 1/A 1/3 as If σ V and σ S depend weakly on A, the 1/A 1/3 dependence makes sense as the leading approximation.A much harder task has been to determine the gluon distribution of nucleons bound in a nucleus, i.e., the nuclear gluon distribution (xg A (x, Q 2 )).The kinematic extension of the the electron -Ion collider (EIC) [21,22] will allow us to examine the non-linear dynamics at low x.When the gluon density becomes sufficiently large at small x, one needs to take into account the effects of gluon recombination (gluon-gluon fusion) leading to nonlinear corrections to the DGLAP evolution equations [23][24][25].Indeed the gluon-gluon recombination processes cause that the growth of the gluon density is slowed down at smaller values of x and Q 2 (but still Q 2 ≫Λ 2 QCD ).In the Gribov-Levin-Ryskin-Mueller-Qiu (GLR-MQ) approach [23,24], the gluon recombination is addressed by analyzing so-called "fan" diagrams, where two gluon ladders merge into a gluon or a quark-antiquark pair.Adding these contributions to the DGLAP equations yields the nonlinear GLR-MQ evolution equations [23,24], where the nonlinear term tames the growth of the PDFs at small x and leads to their suppression.One of the important outcomes of studied in Ref. [26] is the existence of the saturation scale where Q 0 and x 0 are free parameters) which is a characteristic scale at which the parton recombination effects become important.The solution to the non-linear equation has the property of the geometric scaling in the regime where k < Q s (x) whereas in the case when k > Q s (x) the solution enters the linear regime, where k is the gluon transverse momenta.Effects of small-x nonlinear corrections to the DGLAP evolution equations due to gluon recombination have been extensively studied in the literature [27][28][29][30][31][32].Recently in Ref. [33], the authors have considered the nonlinear GLR-MQ evolution equations for nPDFs using the "brute force" method in the momentum space.The authors [33] confirmed the importance of the nonlinear corrections for small x 10 −3 , whose magnitude increases with a decrease of x and an increase of the atomic number A. This paper is organized as follows.In the next section the theoretical formalism is presented, including the GLR-MQ evolution equation.In section 3, we present a detailed analytical analysis and our main results for the nuclear gluon density and predictions of the non-linear effects at higher order accuracy.In the last section we summarize our findings.

Formalism
The nonlinear corrections in the GLR-MQ evolution equations for nPDFs are defined by the following forms and where ∂fi(x,Q 2 ) ∂lnQ 2 | DGLAP for the parton distributions refer to the standard DGLAP evolution equations.Here R is the characteristic radius of the gluon distribution in the hadronic target.R A for a nuclear target with the mass number A is defined by R A = 2 GeV −1 ×A 1/3 [33].The value of 2 GeV −1 depends on how the gluons have a hotspot-like structure within the nucleon.Here χ = x x0 and x 0 is the boundary condition that the gluon distribution joints smoothly onto the linear region.The second terms in the right-hand sides of Eqs.( 7) and ( 8) are expected to become important and related to the recombination of the gluons in the low-x region, when the gluon density is very large.This is known as the phenomenon of gluon saturation.Since the parton distributions in bound and free protons are different, f A (x, Q 2 ) =f (x, Q 2 ), therefore the ratio of structure functions is observed to deviate clearly from unity.The nuclear modifications at x 0.1 are referred to as shadowing.The nuclear structure function F A 2 , in the QCD-improved parton model (in leading order (LO) of α s , or in the DIS-scheme in any higher order), can be written in terms of its parton distributions as where e q is the quark charge, and q A (q A ) is the quark (antiquark) density in the nucleus A. The nuclear structure function, with assumed flavor symmetric antiquark distributions, becomes a summation of valence quark and antiquark distributions The nonlinear equations (i.e., Eqs.( 7) and ( 8)) show that the strong rise that is corresponding to the linear QCD evolution equation at small-x and Q 2 can be tamed by screening effects.After successive integrating of both sides of Eqs.( 7) and ( 8) with respect to ln Q 2 and some rearranging, we find the nonlinear distribution functions in terms of the linear by the following forms and Integrating the first terms in the left and right hands of Eqs.( 11) and ( 12) and using the linear and nonlinear initial conditions xf A i (x, Q 2 0 ) (given by Eqs.( 15), ( 19) and ( 20) below), we find the nonlinear corrections (NLCs) to the parton distribution functions by the following forms 1 For future discussion please see the Appendix.and Here ) are the linear parton distribution functions at the scales of Q2 and Q 2 0 respectively, and obtained from the coupled DGLAP evolution equations using the modified nuclear distribution functions at the initial scale 2 .The initial nuclear parton distributions are provided at a fixed Q 2 (≡Q 2 0 ), due to a free nucleon distribution function, f i (x, Q 2 0 ), and a multiplicative nuclear modification factor, w i (x, A, Z), as The nuclear modification is based on the QCD analysis available in the literature [14,19,[34][35][36][37][38], and assume the following modification function where An advantage of the cubic form with the additional term d i in contrast to a quadratic-type function, i.e., without d i , is that the weight function becomes flexible enough to accommodate both shadowing and anti-shadowing in the valence quark distributions [14].
The nonlinear corrections enter both the gluon and the sea-quark distributions at small x through (i) modifications of the initial distributions and (ii) the presence of additional nonlinear terms in the Q 2 -evolution equations.To study the possible importance of nonlinear corrections, we base our initial gluon and singlet distribution xg N LC (x, Q 2 0 ) and ) by imposing nonlinear corrections on linear distribution functions.The nonlinear corrections to the gluon distribution, at the initial scale Q 2 0 , is obtained from the results in Ref. [39] as where The nonlinear terms in the right-hand side of evolution equations (i.e., Eqs.( 7) and ( 8)) are defined by xg A sat , and this is the value of the gluon which would saturate the unitarity limit in the leading shadowing approximation.In Fig. 1 we show the gluon saturation as a function of the mass number A is expected to occur for various values of Q 2 [40].In Fig. 2, the gluon distribution, α s xg A sat increases as the mass number A increase at the initial scale Q 2 0 = 1.69 GeV 2 [41].Therefore, the effect of the gluon saturation is expected to be larger in heavy nuclei, and to be important for small values of Q 2 .We rewrite Eq.( 17) by using Eq.( 16), to take into account the nonlinear correction to the nuclear gluon distribution at the initial scale for x < x 0 as We note that in Eq.( 19), xg A,NLC →xg A when R A →∞ and xg A sat →∞ and also we see that xg A,NLC →xg A sat when x→0. Moreover xg A,NLC joins smoothly onto xg A at x = x 0 (= 10 −2 ).The nonlinear corrections to the gluon distribution are reflected in the sea-quark distributions q A s (x, Q 2 ) which at small x are predominantly driven by the gluon and modified the nuclear structure function, as 4 The shadowing corrections to the gluon distribution are reflected in the seq-quark distributions which the seq-quark starting distribution in the region x < x 0 in proportion to the shadowing correction to the gluon by the following form [39] xq .
The weight function w i (x, A, Z) for the linear distribution functions can be obtained from the three constrains for the nuclear distributions as the nuclear charge Z, mass number A and momentum conservations5 are defined by the following forms [14,19,[33][34][35][36][37][38]] For a detailed investigation of these functions, we constrain our results to the functions defined in Ref. [35].The gluon distribution at low x is dominant, therefore we used the standard gluon distribution at the input scale Q 2 0 = 1.69 GeV 2 obtained from CT18 set of the free proton PDFs [43], i.e., where the coefficients a 0−4 are listed in Ref. [43].The weight function w g (x, A, Z) for the gluon distribution function is defined by the following form [35] where the coefficients at the next-to-leading order (NLO) and the next-to-next-to-leading order (NNLO) approximations are listed in Ref. [35].The strong coupling is set equal to α s (M z ) = 0.118 for both the NLO and NNLO approximations.In Figs. 3 and 4, we show representations of the nonlinear corrections to the gluon modification functions at the initial scale Q 2 0 = 1.69 GeV 2 for two selected nuclei, C-12 and Pb-208 at the NLO and NNLO approximations, respectively.The nuclear gluon distribution functions are analyzed using the CT18 proton PDF set as a baseline in these figures (i.e., Figs. 3 and 4) [43].The nuclear modification factors have been extracted from QCD fits to the nuclear and neutrino(antineutrino) DIS and Drell-Yan data 6 .The resulting nonlinear corrections to the gluon distribution function are presented in Fig. 5 for carbon (left) and iron (right) at Q 2 0 = 2 GeV 2 in the NLO approximation.To achieve this, we used of the gluon distribution for a free proton defined in Ref. [44] as where the coefficients at the NLO approximation are listed in Refs.[19,44].The weight function for the nuclei of carbon and iron has the same form in Eq. ( 16) which the coefficients are presented in Ref. [19] in which the effects of shadowing, anti-shadowing, fermi motion and the EMC regions are included.In Fig. 6, we compare the nonlinear and linear gluon distributions in lead at the NNLO approximation to those of JR09 [45] at The nuclear gluon distribution is obtained from JR09 parametrization at the input scale by the following form of the free proton PDFs where the parameters in weight function is listed in Ref. [34].To quantify the magnitude of NNLO corrections, we present the nonlinear corrections of nuclear gluon distributions obtained at the input scale of the CT18 and JR09 parametrizations in Figs. of x.Therefore, Eq.( 13) changes to an approximate relation at the NNLO accuracy as In Figs.3-5, we observe that xg A,N LC (x, Q 2 0 ) =xg A (x, Q 2 0 ) at the NLO accuracy, therefore the evolution of the nuclear gluon distribution functions with the nonlinear corrections are defined by the following form For evolution of the nonlinear corrections of the nuclear gluon distributions, we need to a gluon analytical distribution function for a free proton at the scale Q 2 in the NLO and NNLO approximations.In literatures, usually, the gluon analytical distribution function at the LO approximation have been defined.To do it, we extend the analytical solution used in the DGLAP evolution to take into account the nonlinear corrections in Eqs.( 27) and ( 26) at the NLO and NNLO approximations, respectively.We solve the DGLAP evolution equation using Laplace transform techniques in the next section.

Higher order corrections to the gluon distribution
According to the DGLAP Q 2 -evolution equations, the singlet distribution function leads to the following relation of integro-differential equation where P qq and P qg are the quark-quark and quark-gluon splitting functions calculated to the desired order in α s [46][47][48].Here < e 2 > is the average of the charge e 2 for the active quark flavors.Also, < e 2 >= n −1 f n f i=1 e 2 i , and the symbol ⊗ denotes convolution according to the usual prescription.Considering the variable definitions υ≡ ln(1/x) and w≡ ln(1/z), one can rewrite Eq. ( 28) in terms of the convolution integrals and new variables as )) for lead at Q 2 0 = 2 GeV 2 in the NNLO approximation.The delta values are differences between the nonlinear and linear distribution functions where The Laplace transform of H(α s (Q 2 ), υ) , s are given by the following forms Consequently, we can rewrite Eq.( 29) in the Laplace space s, by using the convolution theorem for Laplace transforms and considering the fact that the Laplace transform of the convolution factors are simply the ordinary product of the Laplace transform of the factors, i.e., where and The coefficient functions Φ and Θ in the Laplace space s at the LO approximation are given by Θ (0) where ψ(x) is the digamma function and γ E = 0.5772156... is Euler constant.The explicit expressions for the NLO and NNLO kernels in s space are rather cumbersome; therefore, we recall that we are interested in investigation of the kernels in small x [49][50][51].In the Laplace space, we consider the kernels at small s, as the two and three-loop kernels read (2) with the color factors 3 and T f = 1 2 n f associated with the color group SU (3) and n f being the number of flavors.The strong coupling satisfies the renormalization group equation, which up to NNLO reads where β 0 , β 1 and β 2 are the one, two and three loop correction to the QCD β-function.The standard representation for QCD couplings in NLO and NNLO (within the MS-scheme) approximations have the forms where t = ln Q 2 Λ 2 and Λ is the QCD cut-off parameter [52].Consequently, the discretized form of Eq.( 32) for the gluon distribution reads where the kernels k (ϕ) (α s (Q 2 ), s) and h (ϕ) (α s (Q 2 ), s) contain contributions of the s-space splitting and coefficient functions up the NNLO approximation.These kernels can be evaluated from s-space results by the following forms The inverse Laplace transform of coefficients k(a s (Q 2 ), s) and h(a s (Q 2 ), s) in above equations are defined respectively as kernels and The kernels are dependent on υ and the running coupling at the higher order approximations.In order to obtain an analytical form for these kernels at higher order approximations, we consider the terms of the order 1/s as these terms are dominant at higher order [53].Therefore, we have Consequently, the general analytical expressions for the gluon distribution function in x-space at the higher order approximations are given by Having an analytical proton structure function and its derivative with respect to ln Q 2 , one can extract the gluon distribution function at any desired x and Q 2 values.Using a parameterization suggested by authors in Ref. [54] on the proton structure function in a full accordance with the Froissart predictions [55].The explicit expression for the F 2 parameterization, obtained from a combined fit of the H1 and ZEUS collaborations data [56] in the range of the kinematical variables x and Q 2 ( x < 0.01 and 0.15 < Q 2 < 3000 GeV 2 ), is given by and where Here M and µ 2 are the effective mass a scale factor respectively.The effective parameters in Eq.( 44) are defined in Refs.[54] and [57].

Results for nonlinear nuclear gluon distribution function
Using the analytical approach outlined above (i.e., Eq.( 42)) in the NLO and NNLO approximations, we solve the nonlinear gluon distributions for nuclei at low x as and Now we present our numerical results of the nonlinear gluon distribution for light and heavy nuclei in the x − Q 2 kinematic regions, where the nonlinear corrections are important.The computed results of the nonlinear gluon distribution function for Au-197 compare with the suggested method by Rausch, Guzey and Klasen (the RGK model) [33].This was based on the brute force method, where the authors in Ref. [33] have been extended the numerical algorithm used in the QCDNUM16 DGLAP evolution code [58] to take into account the nonlinear corrections as the nCTEQ15 nPDFs [59] are used as baseline PDFs.
In Fig. 7 we show representations of the nonlinear gluon distribution functions for Au-197 at the scales Q 2 = 4, 16 and 100 GeV 2 as a function of the momentum fraction x to show the effects of the Q 2 evolution.The nuclear weight functions for the gluon are extracted from the suggested method by Khanpour, Soleymaninia, Atashbar Tehrani, Spiesberger and Guzey (the KSASG20 model) [35], where the CT18 nPDFs [43] are used as baseline PDFs.These results are compared to the RGK model [33], where the nCTEQ15 nPDFs [59] are used as baseline PDFs in the nonlinear GLR-MQ evolution equation.In the RGK model, the dashed-dot curves show the results of the upward evolution from Q 2 0 = 4 GeV 2 to Q 2 = 16 and 100 GeV 2 (green curves) and also show the results of the downward evolution from Q 2 0 = 100 GeV 2 to Q 2 = 16 and 4 GeV 2 (purple curves) [33], respectively.The uncertainties, due to the statistical errors of the coefficient functions of the parametrization of the proton structure function [54] and the nuclear modification functions [35], are shown in Fig. 7.For the NLO analysis, the nonlinear nuclear distribution function for the gluon shows an increase as x decreases, which is similar to what one can observe in the analyses by RGK [33].However, the magnitude of these results is slightly differs at different scales, but they are within the uncertainties error bands.As can be seen in the figure, the nonlinear gluon densities come with relatively large error bands at the critical point between the linear and nonlinear (i.e., x = 0.01), reflecting the fact that there are large errors due to the coefficients in the parametrization of the proton structure function.In Fig. 8, the nonlinear gluon distributions for C-12 and Pb-208 at the NLO approximation are considered at Q 2 = 4 and 100 GeV 2 as a function of x as accompanied with their uncertainties.To quantify the magnitude of the nonlinear corrections, we present ratios of nuclear gluon distributions obtained in the nonlinear corrections over those of the linear.Figure 9 quantifies the size of the nonlinear corrections as a function of the mass number A and x for C-12 and Au-197 at Q 2 = 4 and 100 GeV 2 .The difference between the nonlinear and linear evolved gluon densities grows steadily with a decrease of x.This is largest at the smallest values of x and Q 2 and disappears for x = 0.01.The saturation gluon increases as the atomic number increases, therefore the nonlinear/linear ratio decreases as the atomic number increases.As one can see, the nonlinear/linear ratio is slightly larger for light nuclei than for heavy  nuclei, and this effect is, as expected, mainly due to the large gluon saturation values of heavy nuclei.
In Figs.(10) and (11), the nonlinear corrections to the gluon distribution function at the NLO approximation for the nuclei C-12 and Pb-208 are presented at Q 2 = 10 and 100 GeV 2 as a function of the momentum fraction x, respectively.In these figures, our numerical results, which are accompanied with statistical errors, are compared with the linear results based on the KSASG20 (NLO) parametrization [35].The KSASG20 parametrization is a new set of nuclear parton distribution functions (nuclear PDFs) at the NLO and NNLO approximations in perturbative QCD which include the new CT18 PDFs on proton PDFs.As can be seen in these figures, the effects of nonlinear corrections are noticeable at small x values and the strong growth of gluon distributions are tamed by shadowing effects as x decreases.The solid curves represent the effect of shadowing correction for R A = 2A 1/3 GeV −1 presented  by using Eq.( 45).As can be observed, the nuclear gluon distributions increase as x decreases, which corresponds with the perturbative QCD fits at small x, but these behaviors are tamed with respect to nonlinear terms at the GLR-MQ equation.These tamed behaviors of nuclear gluon distributions due to the shadowing corrections satisfy the Froissart bound in the perturbative QCD means.Hence, as one can see from Figs.10 and 11, deviations from the linear nuclear gluon distributions based on the KSASG20 (NLO) parametrization increases as x decreases.The deviations from the KSASG20 (NLO) nPDFs increases as Q 2 increases and decreases as atomic number A increases (indeed the nonlinear nuclear gluon distributions increases as atomic number A increases), and significant effects are found for heavier nuclei, such as lead.These behaviors for the nonlinear nuclear gluon distributions are similar to the analysis of RGK [33].

Summary
In conclusion, we have studied the effects of adding the nonlinear corrections to the gluon distribution function for light and heavy nuclei at small x analytically.We used the parametrization of the proton structure function to take into account an analytical solution for the gluon density at low x in the NLO approximation.The nuclear modification factors are obtained from KSASG20 nuclear PDFS which are based on the CT18 framework.The shadowing effects of the gluon distribution at small x through modifications of the starting distributions and the presence of additional nonlinear terms in the initial point Q 2 0 at the NLO and NNLO approximations for light and heavy nuclei considered.We obtained the nonlinear corrections for small x in a wide range of Q 2 values.These results show that the magnitude of the nonlinear corrections increases with a decrease of x and an increase of the atomic number A. Our results are consistent, within uncertainties, with the determination of nuclear gluon distribution with the upward and downward evolution from the RGK model, which is based on the nCTEQ15 nPDFs as input.Our determination of nuclear gluon distributions includes error estimates obtained with respect to the coefficient errors in the parametrization of the proton structure function and the nuclear modification function errors.We found differences between our and RGK results at the NLO accuracy, which occur in different assumptions such as the input prameterizations and also the approximate relation between the gluon distribution and the proton structure function due to the Laplace transform method at low x.These results for the nonlinear corrections to the nuclear gluon distribution function may be important for future experiments at the Electron-Ion Collider [21,22], LHeC Collaboration or a Future Circular Collider (FCC) study group [5] and Electron-Ion Collider in China (EiCC) [60] at low x.
where the higher dimensional gluon term G HT is here assumed to be zero.The standard DGLAP evolution equation for singlet and gluon distributions has the following forms: ∂xg(x, Q where P , ij s are the splitting functions in the desired order in α s .

FIG. 9 :
FIG. 9:The ratio of the nonlinear/linear gluon distributions for C-12 and Au-197 at Q 2 = 4 and 100 GeV 2 as a function of x with their uncertainties.