Studies on gluon evolution and geometrical scaling in kinematic constrained unitarized BFKL equation: application to high-precision HERA DIS data

We suggest a modified form of a unitarized BFKL equation imposing the so-called kinematic constraint on the gluon evolution in multi-Regge kinematics. The underlying nonlinear effects on the gluon evolution are investigated by solving the unitarized BFKL equation analytically. We obtain an equation of the critical boundary between dilute and dense partonic system, following a new differential geometric approach and sketch a phenomenological insight on geometrical scaling. Later we illustrate the phenomenological implication of our solution for unintegrated gluon distribution $f(x,k_T^2)$ towards exploring high precision HERA DIS data by theoretical prediction of proton structure functions ($F_2$ and $F_L$) as well as double differential reduced cross-section $(\sigma_r)$. The validity of our theory in the low $Q^2$ transition region is established by studying virtual photon-proton cross-section in light of HERA data. We also investigate the claim that the longitudinal structure function $F_L$ is a sensitive probe of the gluon distribution and extract colinear gluon density $xg(x, Q^2)$ from $F_L$. Finally, we portray a dipole picture of kinematic constraint improved unitarized BFKL equation and examined $r$ (transverse size of the dipole) dependence of dipole cross-section.


Introduction
In perturbative QCD the two well-known parton evolution equations DGLAP [1][2][3] and BFKL [4,5] serve as the basic tools for prediction of parton distribution functions (PDFs) on their respective kinematic domains. The DGLAP approach is valid for the large interaction scale Q 2 since it sums higher order α s contributions enhanced by the leading logarithmic powers of ln Q 2 . On the other hand, BFKL evolution deals with the small-x and semi-hard Q 2 kinematic region by summing the leading logarithmic contributions of (α s ln 1/x) n . However, both the linear evolutions turn out to be inadequate to comply with the unitary bound or the conventional Froissart bound [6]. To restore the unitarity, a higher order pQCD correction is required to shadow the rapid parton growth which eventually leads to the saturation phenomena at very small-x. In this respect, several nonlinear equations have been proposed in recent years to entertain this shadowing correction in DGLAP and BFKL equation.
In pQCD interactions, the origin of the shadowing correction is primarily considered as the gluon recombination (g + g → g) which is simply the inverse process of gluon splitting (g → g + g). The modification of gluon recombination to the original DGLAP equation was first performed in GLR-MQ equation proposed by Gribov-Levin-Ryskin and Mueller-Qui [7][8][9]. It is considered that the interferant cut diagrams of the recombination amplitude are the origin for the negative shadowing effects in the gluon recombination processes [10,11]. The GLR-MQ equation was originally derived using AGK-cutting rules [12] to compute the contributions from these interference processes. However, GLR-MQ does not possess any antishadowing correction term since in this equation momentum conservation is violated in the gluon recombination processes. The shadowing and antishadowing effects are expected to happen in a complementary fashion which in fact ensure the momentum conservation during the process of gluon recombination [13]. In consequence of gluon recombination, a part of gluon momentum lost due to shadowing is believed to be compensated in terms of new gluons with comparatively larger x. This causes an enhanced density of larger momentum gluons responsible for the antishadowing effect. In this respect, to take care of momentum conservation and antishadowing effects in gluon evolution, a modified DGLAP equation was proposed in [11] to replace the GLR-MQ equation. In the derivation of the modified DGLAP equation, the author has used TOPT-cutting rules based on time ordered perturbation theory (TOPT) to calculate the contribution of the gluon recombination, instead of using AGK-cutting rules unlike in the derivation of GLR-MQ equation. This equation naturally comes with two separate nonlinear terms, quadratic in gluon distribution function which are identified as the positive antishadowing and negative shadowing components in the gluon evolution. This suggests that negative shadowing effects are suppressed by antishadowing effects in the gluon evolution process.
However, above mentioned modified DGLAP approach cannot predict the evolution of unintegrated gluon distribution because DGLAP is based on colinear factorization scheme while the later deals with k T -factorization. In this respect, a new unitarized BFKL equation was proposed by Ruan, Shen, Yang and Zhu [14] incorporating both shadowing and antishadowing correction. This modified BFKL (MD-BFKL) equation allows one to study antishadowing effects in unintegrated gluon distribution platform.
However, there are several other BFKL based nonlinear evolution equations that consider corrections from gluon fusion. One of the most widely studied models is Balitsky-Kovchegov (BK) equation [15,16] which is originally derived in terms of scattering amplitude. It can be explained by the dipole model in which the nonlinear term comes from the dipole splitting while the screening effect origins from the double scattering of the dipole on the target. The solution of the BK equation suggests an intense shadowing leading to the so-called saturation of the scattering amplitude showing a complete flat spectrum. However, there is no clear relationship between scattering amplitude and unintegrated gluon distribution so far and this indeed turns out to be model-dependent. One such modeldependent transformation of BK equation for the unintegrated gluon distribution can be found in [17]. Study of the analytical solution of this transformed equation is difficult due to the complicated form of its nonlinear term. Moreover, like GLR-MQ equation the BK equation is also irrelevant to the gluon antishadowing. On the other hand beside inclusion of antishadowing correction term, the nonlinear terms in MD-BFKL equation hold a comparatively simple form which is quadratic in unintegrated gluon distribution as well as running coupling constant with some constant coefficients. These unique features of MD-BFKL equation motivate our current studies on small-x physics.
A numerical solution of MD-BFKL equation suggests a sizable impact of antishadowing effect on gluon evolution particularly in the pre-asymptotic regime [14]. In the literature [18] the phenomenology of the MD-BFKl equation is extended to study the implicit nuclear shadowing and antishadowing effects. In our current work, we try to construct a modified evolution equation by implanting the so-called kinematic constraint or consistency constraint on the MD-BFKL equation. The kinematic constraint (see Fig. 1 ) is implemented in different forms: LDC [19,20], (1.1) [21]. (1.3) This constraint arises as a consequence of BFKL multi-Regge kinematics which suggests the exchanged gluon virtuality is dominated by transverse components while the longitudinal components of the gluon momentum are required to be small i.e. k ′ 2 ≈ k ′ 2 T . The kinematic constraint gives an implicit cutoff on k ′ 2 T as depicted by (1.3). The inequality (1.3) can be considered as a special case of (1.1) recalling the fact that for a given value of k 2 T , a high q 2 T implies an equally high k ′ 2 T [21]. Although there are other constraints coming from energy-momentum conservation [22], this bound is considerably tighter than the later [21]. Besides the introduction of this cutoff on the upper limit of integration found to preserve the scale invariance of the BFKL equation.
Our primary goal of this work is to investigate the nonlinear effects on gluon evolution in terms of MD-BFKL equation supplemented by the kinematic constraint. Our studies are particularly focused on the near saturation region where shadowing effects are dominant over antishadowing effects. The content of the paper is organized as follows. In Sect. 2 we present our kinematic constraint improved MD-BFKL equation starting from the construction of the equation and a brief discussion on some of its main features (Sect. 2.1). We suggest an analytical solution to this equation and sketch a particular solution in terms of x and k 2 T evolution followed by three-dimensional realization of gluon evolution in x-k 2 T phase space (Sect. 2.2). Our theoretical prediction of unintegrated gluon distribution is compared with TMDlib [23] data sets viz. PB-NLO-HERAI+II-2018 [24,25] and ccfmset-A0 [26]. In Sect. 3 we suggest a differential geometric approach towards finding the equation of the critical boundary between high and low gluon density regime. In this section, a phenomenological insight of the geometrical scaling is discussed as well.
In Sect. 4 we explore the phenomenological implication of our solution for unintegrated gluon distribution towards prediction of DIS structure function F 2 and longitudinal structure function F L at HERA. We show the comparisons of reduced cross-sections (Sect. 4.1) as well as virtual photon-proton cross-sections in the transition region (Sect. 4.2) with HERA H1 and ZEUS data. We then investigate the claim that the longitudinal structure function F L is a sensitive probe of the gluon distribution and perform the extraction of collinear gluon density, xg(x, Q 2 ) from our theoretical prediction of F L (Sect. 4.3). The results are compared with HERA as well as various LHAPDF global data-fits viz. NNPDF3.1sx [27], HERAPDF 15 [28], ABMP [29], ATLAS [30] and CT14 [31]. In Sect. 5 we recast our KC improved MD-BFKL equation in terms of dipole scattering amplitude. In this dipole picture, we show that the ordinary MD-BFKL equation reduces to the BK equation while KC improved BFKL equation reduces to the KC improved BK equation [32]. We also investigate r (transverse size of the dipole) dependence of the dipole cross-section and compared our result with that of well-known GBW parametrization. Finally in Sect. 6 we summarize and outline our conclusion as well as possible future prospects.
2 Kinematic constraint improved MD-BFKL equation and its nonlinear effects

Construction of kinematic constraint improved MD-BFKL
The modified BFKL (MD-BFKL) equation reads [14], where f x, k 2 T denotes the gluon distribution unintegrated over the transverse momentum of gluon k T and k 2 0 is the infrared cutoff of the evolution. The first part of (2.1) linear in is the BFKL kernel at leading logarithmic of 1/x (LLx) accuracy. A diagrammatic representation for probing gluon content of the proton at high photon virtuality Q 2 is sketched in Fig. 1 (a). The BFKL kernel corresponds to the sum of gluon ladder diagram Fig. 1(b) generated by squaring the amplitude of Fig. 1(a). The 1st term in the BFKL kernel, involving f (x, k ′ 2 T ) corresponds to diagrams with real gluon emission while the second term takes care of diagrams with virtual corrections. Note that the apparent singularity that is observed at k ′ 2 T = k 2 T cancels between real and virtual contributions. The quadratic terms in (2.1) viz.
c −1 f 2 (x, k 2 T ) depict shadowing and antishadowing correction to the original BFKL equation respectively. The factor πR 2 represents the transverse area populated by gluons. The radius parameter R arises in the QCD cut diagram coupling 4 gluons to 2 gluons (see Fig. 2), the value of which depends on how exactly the gluon ladders couple to hadron. If the gluon ladder couples to two different constituents of hadron, R is characterized by hadronic radius i.e. R ≈ R H = 5 GeV −1 whereas if we consider the possibility of ladder coupling to the same parton then the appropriate choice of R is the radius of valence quark i.e. R = 2 GeV −1 . In the latter case (R = 2 GeV −1 ) we have "hot spots" [8,33].
The interpretation of the real emission term in the BFKL kernel is that we have a splitting k ′ → k + q inside hadron resulting an infinite chain of reggeized gluons labeled as k 1 , k 2 , k 3 , ...k n ( Fig. 1(a)). In high energy limit, the longitudinal components of the gluon momentum are strongly ordered while there is no ordering on the transverse components of the gluon momentum i.e.
In addition, the small-x regime where the BFKL is valid, the gluon virtuality along the chain must be dominated by the transverse components of the gluon momentum, The above kinematics corresponding to (2.2a),(2.2b) and (2.3) is referred as multi-Regge kinematics. The inequality (2.2a) implies z = x x z = x n−1 xn ≪ 1 and since transverse momenta are of the same order: which is so-called kinematic constraint or consistency constraint for real gluon emission [17,19,21,34]. The BFKL kernel in (2.1) is at leading logarithmic 1/x (LLx) accuracy. However higher order corrections to the BFKL equation are already been evaluated up to NNLx accuracy, which turns out to be quite large from the analysis in [35]. Implementation of the constraint (2.4) on the evolution of BFKL based equations makes the theory more realistic in the sense that this ensures the participation of only the LLx part of higher order correction in the evolution. This, in fact, portrays the importance NNLx correction in BFKL kernel.
Recalling that BFKL equation can be written as an integral equation of f (x, k 2 T ) [7,36], the kinematic constraint (2.4) can be imposed onto the real emission part of the BFKL (2.8) In derivation of (2.8) we have neglected the term x ∂f (0) ∂x . This is justified in the sense that x ∂f (0) ∂x is much less singular than x ∂f ∂x at small-x. Moreover if we consider f (0) is particularly contributed by the diagrams shown in Fig. 1 [19,36]. To simplify the distribution function f ( T ) corresponding to real emission term in (2.8) we have incorporated Regge like behavior of gluon distribution. Before the advent of pQCD, Regge theory is been considered as a very successful theory regarding the phenomenological analysis of total hadronic and photoproduction cross-sections. These non-perturbative processes can economically described by two Regge contributions namely pomerons with intercept slightly above unity (α P (0) ≈ 1.08) and Reggons with intercept α R (0) ≈ 0.5. The Regge behavior corresponding to sea quark and anti-quark distributions is given by q sea ∼ x −α P whereas that of valence quark distribution is given by q v ∼ x −α R . On the other hand in pQCD also, particularly the small-x region is considered to have higher possibility towards exploring Regge limit of pQCD. Note that the BFKL dynamics itself is based on the concept of pomeranchuk theorem or pomeron: the Regge-pole carrying the quantum-numbers of the vacuum. However, the BFKL pomeron (also called hard pomeron) should be contrasted with non perturbative description of pomeron (soft pomeron) in the sense that they pose relatively different magnitude of intercept. For BFKL-pomeron the intercept is where λ BFKL = 3αs π 28ζ(3), ζ being Reimann zeta function [37]. The typical value of λ BFKL for α s =0.2 is ∼ 0.5 implies α BF KL P (0) ≈ 1.5 which is potentially large in magnitude compared to soft pomeron (α P (0) = 1.08).
In pQCD the small-x behavior of the parton distribution is supposed to be controlled by intercept of appropriate Regge trajectory. Regge model provides parametrizations of DIS distribution functions, f i (x, Q 2 ) = A i (Q 2 )x −λ i (i= (singlet structure function)and g (gluon distribution)), whereλ i is the pomeron intercept minus one (α P (0) − 1). At small-x, the leading order calculations in ln(1/x) with fixed α s predicts a steep power law behavior of f x, k 2 ∼ x −λ BFKL [22,38]. This motivates us to consider a simple form of Regge factorization as follows, (2.11) The above equation (2.11) is our kinematic constraint improved MD-BFKL equation in LLx accuracy. Before going into in depth phenomenological analysis and consistency of the equation towards experimental result, we want to highlight two important features of the evolution equation (2.11). As we follow, in the region far below saturation limit, for the canonical choice of λ ∼ 0.5, the quadratic term in (2.11) tends to vanish since (1 − 2 2λ−1 ) → 0. This indicates that below saturation region both the contribution from shadowing correction ∂f shad (x,k 2 T ) ∂ ln 1 x and antishadowing correction ∂f antishad (x,k 2 T ) ∂ ln 1 x coexists but seem to balance each other for which nonlinear effect becomes negligible. Therefore in the below saturation regime, we can revert (2.11) back to the original BFKL equation choosing λ = 0.5.
On the other hand, in the vicinity of saturation limit, our interpretation of gluon distribution for the nonlinear term has to be reviewed. In this region, the gluon distribution becomes flat which makes our Regge factorization for the nonlinear term in (2.10) invalid. Rather essentially we should take the approximation which is supposed to justify our understanding of the saturation region. Taking this approximation in (2.8) one can arrive at (2.12) Equation (2.12) serves the evolution of gluon near the saturation limit. Here the λ dependence of the nonlinear term has been completely wiped up, unlike (2.11). Note that in this regime, the shadowing contribution becomes twice as that of anishadowing effect forecasting a net shadowing effect in the evolution. Another important features of (2.12) is that it can be reduced to Balitsky-Kovchegov (BK) equation [15,16] (with KC) by simple scale translation prescription which will be demonstrated in the latter part of this text.

Analytical solution of KC improved MD-BFKL
In this section, we present an analytical solution of our KC improved MD-BFKL equation, particularly in the saturation region where shadowing effect is the dominant one. Since our calculations are limited to fixed strong coupling α s , so do fix λ, therefore the solution for both (2.11) and (2.12) will exhibit the same form only differ by the coefficients of their respective quadratic terms.
Recalling that there is no ordering for transverse momenta k ′ T ⋍ k T in BFKL multi-Regge kinematics, this allows us to write the gluon distribution in Taylor series, T denotes the higher order terms of the series. Above series is expected to be well-behaved as long as no ordering condition, k ′ T ⋍ k T is concerned and the higher order terms O become insignificant. Now neglecting these higher order terms we can express (2.13) as T is a standard hypergeometric function. In the above calculations we have taken infrared cutoff k ′ 2 T min = 1 GeV 2 since for unified BFKL-DGLAP framework this provided a very consistent result towards HERA DIS data for proton structure function F 2 [45]. In (2.23) (k 2 T ) λ /λ is the only dominant term since other terms are significantly small. The contribution from the logarithmic term in (2.23) is negligible in comparison to the net contribution for all k 2 T . However phenomenological studies shows the term involving hypergeometric function in (2.23) becomes irrelevant towards change in k 2 T i.e. it posses a constant value (≈ −4.79) (see Fig. 3(a)). This constant contribution is insignificant to the net contribution for ζ(k 2 T ) if k 2 T is enough high. But for small k 2 T this contribution can not be neglected since at this range, the k 2λ T /λ contribution itself is small. In consequence, this constant contribution can be treated as a small perturbation ǫ to the dominant term k 2λ T /λ . We have performed phenomenological determination of the constant perturbation parameter ǫ using standard non-linear regression method (see Fig. 3(b)). Now we are set to solve our original equation (2.14) which is indeed a 1st order semilinear (nonlinear) PDE. Our analytical approach of solving the same involves two steps: first we will express the nonlinear PDE in terms of a linear PDE then we will solve the linear PDE via. change of coordinates.
which is in fact a linear PDE in ∂ω ∂x , ∂ω and ω. Now we construct a new set of co-ordinate σ ≡ σ(x, k 2 T ) and η ≡ η(x, k 2 T ) such that it transforms (2.24) into an ODE. To be specific we shall define this transformation in such a way that it is one to one for all (x,y) in some set of points D in xy plane. This will allow us to solve (2.24) for x and y as functions of σ and η. The only requirement is that we should ensure the Jacobian of the transformation Next we want to recast (2.24) in (σ, η) plane computing the derivatives via chain rule: Now since we want above equation to become an ODE, we require either, If we consider xσ x + ξσ k 2 T = 0, the solution σ is constant along the curves that satisfy where l = − αsNc πλ and C 1 is the constant of integration. Now we can choose the new We rewrite (2.25) as is a standard Gamma function and G e − k −2λ lλ x is an arbitrary continuously differentiable function. The parameters m = ǫαsNc π and n = αsNc πλ are coming from (2.23). Now recalling that f = ω −1 we get the unintegrated gluon distribution , (2.30) which is the general solution of the KC improved MD-BFKL equation. In the following section, we try to get particular solutions for the PDE applying some initial boundary condition and present an analysis of the phenomenological aspects of (2.30).

x and k 2
T evolution In this section, we study the small-x dependence of gluon distribution by picking an appropriate input distribution at some high x, as well as the k 2 T dependence of gluon distribution setting input distribution at some low k 2 T . The input distributions are the non-perturbative inputs, that perturbative QCD cannot predict.
Let us rewrite (2.30) rearranging a bit T lλ x is an arbitrary continuously differentiable function. We will try to evaluate the functional form of this function applying some initial boundary condition on it. For k 2 T evolution we set the initial distribution at (x, where F 0 is the unintegrated gluon distribution at (x, k 2 0 ). Equation (2.32) is the functional form of G. Replacing t with any other argument will give us the value of G at that particular argument. We put t = e − k −2λ T lλ x in (2.32) which gives us the l.h.s. of (2.31) and then we solve (2.31) for f (x, k 2 T ). The solution for k 2 T evolution turns out to be Similarly setting the initial distribution at (x 0 , k 2 T ), we obtain the solution for x-evolution as follows We have plotted the solution for both x evolution (2.33) and k 2 T evolution (2.34) and compared our result with TMDlib [23] Datasets viz. PB-NLO-HERAI+II-2018 [24,25] and ccfm-set-A0 [26] in Fig. 4-5. The input gluon distribution F 0 in our calculation is taken from PB-NLO-HERAI+II-2018. The PB-NLO-HERAI+II-2018 parameterization includes collinear and TMD parton densities obtained from fits to the HERA precision cross-section measurement. The PDFs are evolved by NLO DGLAP evolution using a recently developed theory called Parton Branching (PB) method [24,46] which allows the calculation of both collinear and TMD densities simultaneously. On the other hand the other parameterization ccfm-set-A0 [26] includes uPDF sets extracted from CCFM evolution taking starting distribution from GBW parametrization [47]. Note that our comparison with Data is a little indirect one since the TMDlib parameterizations provides datasets in terms of dipole scattering amplitude φ(k T , x). Although the role of dipole scattering amplitude φ(k T , x) is same as the unintegrated gluon density f (x, k 2 T ) in an evolution equation, they are differ by definition. However no clear relationship between φ(k T , x) and f (x, k 2 T ) is found so far, rather it turns out to be model dependent [48][49][50]. It is expected that φ(k T , x) and f (x, k 2 T ) differ only by some simple scale translation f ∼ k 2 T φ [14,[48][49][50]. In our case we take a more convenient relation between f (x, k 2 T ) and φ(k T , x) as follows The phenomenology behind the relation (2.35) is discussed in Sect. 5. Here φ(k T , x) preserves its dimensionless behavior since dimension of R 2 Another important point is that the gluon density provided by the TMDlib are hard scale µ 2 dependent while our formalism is µ 2 scale independent. However in our kinematic region of moderate k 2 T it is justified to compare with those data sets if we take care of the BFKL kinematics k 2 T ≈ µ 2 while extracting data. This may seem to overestimate the comparison by magnitude but satisfactory information about the growth of gluon distribution can be expected towards small-x and moderate k 2 T . In Fig. 4-5 we have shown k 2 T and x evolution obtained for two different form of shadowing: conventional R = 5 GeV −1 (order of proton radius) where gluons are spread throughout the nucleus and extreme R = 2 GeV −1 where gluons are expected to concentrated in hotspots within the proton disk recalling that πR 2 is the transverse area within which gluons are concentrated inside proton. From the Fig. 4-5 it is clear that shadowing correction are more prominent when gluons are concentrated in hotspots within proton.
The k 2 T evolution in Fig. 4 is studied for the kinematic range 1 GeV 2 ≤ k 2 T ≤ 100 GeV 2 corresponding to six different values of x as indicated in the figure. The input distribution F 0 (x, k 2 0 ) is taken at k 2 0 = 1 GeV 2 from PB-NLO-HERAI+II-2018 dataset. Our evolution shows a similar growth as datasets for all x at conventional shadowing condition R = 5 GeV −1 . The extreme shadowing form R = 2 GeV −1 seems to overestimate the dataset particularly at very small-x (≤ 10 −4 ) values. It is also observed that the growth of f (x, k 2 T ) is almost linear for the entire kinematic range of k 2 T . This is expected since the net shadowing term in (2.12) has 1/k 2 T dependence, which suppresses the contribution from the shadowing term at large k 2 T . The x evolution of unintegrated gluon distribution f (x, k 2 T ) is shown in Fig. 5 for six different values of k 2 T ranging from 4 GeV 2 to 100 GeV 2 . The input is taken at higher x value x = 10 −2 and then evolved down to smaller x value upto x = 10 −6 thereby setting the kinematic range of evolution 10 −6 ≤ x ≤ 10 −2 . We observed that the singular x −λ growth of the gluon is eventually suppressed by the net shadowing effect. However it is hard to establish the existence of shadowing for x ≥ 10 −3 . The obvious distinction between the two form of shadowing R = 5 GeV −1 (conventional) and R = 2 GeV −1 ("hotspot") is also observed towards small-x (x ≤ 10 −3 ). However, at large k 2 T values (k 2 T = 36 GeV 2 , 64 GeV 2 and 100 GeV 2 ) the distinction between the two form of shadowing becomes very small even at x = 10 −6 . Interestingly towards very small-x (x ≤ 10 −5 ), at small k 2 T values (k 2 T = 4 GeV 2 , 16 GeV 2 and 25 GeV 2 ) the gluon distribution becomes almost irrelevant of change in x. This could be a strong hint for the possible saturation phenomena in the small-x high density regime.
In the two different forms of the MD-BFKL equation from linear BFKL equation reflects the underlying shadowing correction. It is also observed that shadowing effect is more intense in MD-BFKL with KC than MD-BFKL without KC.

Complete solution of KC improved MD-BFKL
In this section we implant a functional form of the input distribution (more likely a dynamic one) on the general solution of our KC improved MD-BFKL equation (2.30) and try to obtain a parametric form of the solution. The underlying motivation towards doing so is that this allows us to evolve our solution for both x and k 2 T simultaneously in x-k 2 T phase space which help us to portray a three-dimensional realization of the gluon evolution.
Recall the well-known solution of the linear BFKL equation [37] f (x, where λ = 3αs π 28ζ(3), ζ being Reimann zeta function and Ω = 32.1α s . The nonperturbative parameter k 2 s = 1 GeV 2 and the normalisation constant β ∼ 0.01 [14]. Since far below saturation region both linear and nonlinear equation should give the same solution, therefore we can take (2.36) as the input distribution for our general solution of KC improved MD-BFKL equation.
First we try to find the functional form of the arbitrary differentiable function G e − k −2λ (2.37) Setting initial parameters at (x 0 , k 2 T ) we get initial distribution (2.36) as We denote the argument of G at (x 0 , k 2 where, k 2 T = (−lλ ln(τ x 0 )) −1/λ . Note that (2.39) is the functional form of G. We substitute T lλ x in (2.39) which gives us the l.h.s. of (2.37) and then solve (2.37) for f (x, k 2 T ), where, For simplicity let us denote r.h.s. of (2.40) by γ i.e. .

(2.43)
From (2.36) we have the input distribution   .45) itself, therefore, we do not have to depend on the experimental data fits for input distribution unlike we did in the previous section of x and k 2 T evolution. This parametric form of the solution actually helps us to explore the three-dimensional insight of the gluon evolution in x-k 2 T phase space. However using (2.45) one may also study x and k 2 T evolution separately by setting x = x 0 for k 2 T evolution or k 2 T = k 2 0 for x evolution. In Fig. 7(a-b) we have shown a comparison between TMDlib [23] data set PB-NLO-HERAI+II-2018 [24,25] (Fig. 7(a)) and our solution of KC improved MD-BFKL equation (Fig. 7(b)) in three dimension. The kinematic region for our study is set to be 10 −6 ≤ x ≤ 10 −2 and 1 GeV 2 ≤ k 2 T ≤ 10 3 GeV 2 . From Fig. 7 it is clear that our result exhibits similar behavior as of the dataset. However, a small disagreement between the two surfaces of the datasets and theory is observed towards very high k 2 T particularly for the range 200 GeV 2 ≤ k 2 T ≤ 10 3 GeV 2 . In that range of k 2 T our theoretical prediction shows a linear rise in gluon distribution with no significant shadowing suppression. This is expected because, as mentioned earlier in k 2 T evolution, the shadowing term is proportional to the factor 1/k 2 T , which suppresses the net shadowing effect at larger values of k 2 T . Note that even in that range 200 GeV 2 ≤ k 2 T ≤ 10 3 GeV 2 at very small-x i.e. 10 −6 ≤ x ≤ 10 −4 our result seems to regain the compatibility towards the dataset.
In Fig. 8 we have shown density plots of f (x, k 2 T ) in x-k 2 T domain to examine the sensitivity of f (x, k 2 T ) towards the parameter λ. Density plot allows us to visualize and distinguish the kinematic regions with high/low f (x, k 2 T ) in x-k 2 T plane which turns out to be more informative then any ordinary 3D plot. In Fig. 8 the plots are sketched for two canonical choices of λ viz. λ = 0.4 GeV −1 and 0.6 GeV −1 corresponding to two α s values 0.15 and 0.23. Our solution seems to be very sensitive towards a small change in λ. An apparent 20% change in λ (0.4 to 0.6) suggests approximately around one order of magnitude rise in gluon distribution f (x, k 2 T ) for an approximate limit of x and k 2 T : 10 −6 ≤ x ≤ 10 −5 and 50 GeV 2 ≤ k 2 T ≤ 100 GeV 2 . It is also observed that the range of high k 2 T and very small-x is the high gluon distribution f (x, k 2 T ) region where gluons are mostly populated. In Fig. 9 we have shown R sensitivity of our solution for two choices of the shadowing parameter R viz. R = 5 GeV −1 and 2 GeV −1 . A satisfactory shadowing effect is observed from the comparison of the two plots. The extreme form of shadowing (R = 2 GeV −1 ) is found to suppress atleast 10-20% magnitude of gluon density than the conventional shadowing (R = 5 GeV −1 ) in the high k 2 T and small-x region.

Equation of Critical line and prediction of saturation scale: a differential geometric approach
The so-called saturation physics allows one to study the high parton density region in the small coupling regime. The transition from the linear region to saturation region is characterized by the saturation scale. The saturation momentum scale Q s is the threshold transverse momentum for which non-linearity becomes visible. The boundary in (x, Q 2 ) or (x, k 2 T ) plane along which saturation sets in is characterized by the critical line. An important feature of our analytical solution to the KC improved MD-BFKL equation is the finding of the equation of the critical line which is supposed to mark the boundary between dilute and dense partonic system in x-k 2 T phase space. Although the gluon saturation can be achieved only when Q s ∼ Q (or k 2 T ), the observables already become sensitive to Q s during the approach to saturation regime. This property is known as geometrical scaling in DIS inclusive event which means that instead of depending on k 2 T and x separately, the gluon distribution depends on a single dimensionless variable In recent years, this geometrical scaling property of DIS observables is studied very extensively for various frameworks [32,[51][52][53][54]. In [51] an analysis of the saturation scale has been performed in the platform of resummed NLLx BFKL where the saturation scale was calculated via the relation which has been repeatedly derived in several literature [10,55,56]. In [32] the saturation scale Q s was obtained from the numerical solution of a nonlinear equation by finding the maximum of the momentum distribution of the gluon. Another approach for determination of Q s can be found in [17] where a parameter β is defined as the relative difference between the solutions to the linear and nonlinear equation, .

(3.3)
The crucial parameter β actually depicts the percentage deviation that the non linear solution shows from the linear one and it lies in the order of 0.1-0.5 (or 10-50% deviation). Our approach towards studying geometrical scaling and critical line is primarily based on the basic understanding from differential geometry in particular gradient of a function which is considered as the direction of steepest ascent of that function. Each component of the gradient gives us the rate of change of the function with respect to some standard basis i.e. it gives us an idea about how fast our function grows or decays or saturates with respect to the change of the variables. One important advantage of choosing gradient is that it is a two dimensional object since it does not posses any component along f (x, k 2 T ) axis in R 3 . This ensures that it does not have any direct dependence on the magnitude of gluon density f (x, k 2 T ), rather it depends on the rate at which f (x, k 2 T ) changes with respect to x and k 2 T change. This property of gradient actually helps us in distinguishing out the saturation region and linear region although the distribution function f (x, k 2 T ) has large variation in order of magnitude for different regions in x-k 2 T phase space. Recall that towards small-x, gluon evolution is suppressed due to shadowing effect as sketch in Fig. 7 which motivates us to evaluate the gradient of f (x, k 2 T ) particularly along x basis. For simplicity we consider an unit vector ν alongη (= 1/x) basis, we can project the gradient ∇f (η, k 2 T ) along ν via dot product ∇f (η, k 2 T ). ν. This scalar quantity ∇f (η, k 2 T ). ν can also be interpreted as the directional derivative along the direction ν, Taking the Euclidean norm yields where the negative (-) sign is for the descending function (or negative slope). We obtained a family of contours (or level curves)η(k 2 T ) inη-k 2 T plane solving (3.5) as shown in Fig. 10(a). Each contour depicts a constant gradient g along the curves. The set of contours can also be identified as some set of possible saturation scales. As sketch in Fig. 10(a) theη-k 2 T plane is divided into two regions: low gluon density region where the spacing between two consecutive contours is very small and high gluon density region where the spacing becomes very large compared to previous. This distinction in contour spacing  [47].
between the two regions comes from the fact that the gradient changes very fast until the saturation is reached and then after reaching the saturation boundary gradient changes very slowly or almost freezes for further increase inη(k 2 T ) as can be seen in Fig. 10(a). In high density region the contour curves tend toη(k 2 T ) → ∞ towards high k 2 T . For low gluon density region, shadowing effects are negligible and the contours become almost parallel straight lines.
Let us try to find the equation of critical line which divides the two regions inη-k 2 T space. The level curves of the function g = ±|∇ ν f (η, k 2 T )| are two dimensional curves that can be obtained by setting g = k where k is a constant (k ∈ R). Therefore the equations of the level curves are given by Now for a known initial saturation momentum scale Q s0 (1/x 0 ) we can predict equation of the critical line using (3.5) and (3.6). The equation of the critical line is the equation for the level curve of the function g = ±|∇ ν f (η, k 2 T )| that passes through the point (Q s0 (η 0 ),η 0 ). First we find the value of k by plugging the point (Q s0 (η 0 ),η 0 ) into (3.6) Now the level curve passing through (Q s0 (η 0 ),η 0 ) is obtained by setting which is the equation of the critical boundary. The knowledge of an appropriate initial saturation scale Q s0 (1/x 0 ) allows one to separate out the linear and saturation region using (3.8). In Fig. 10(b) we have sketched a possible critical line obtained from (3.8) for the choice of initial saturation scale Q 2 s0 (η 0 ) ⋍ 2.8 GeV 2 at η 0 = 10 6 (or x 0 = 10 −6 ) which is taken from the calculation from the original saturation model by Golec-Biernat and Wusthoff [17,47]. A rough agreement between our prediction and that of GBW model is observed in Fig. 10. However Q s given by (3.8) is found to have weaker x dependence than the one from GBW model Q ′ 2 s . The saturation scale has direct dependence on partons per unit transverse area. Smaller x suggests larger parton density giving rise to a larger saturation momentum scale, Q 2 s . In other words the saturation scale Q s depends on x in such a way that with decreasing x one has to probe smaller distances or higher Q 2 in order to resolve the dense partonic structure of the proton which is clear from our analysis.

DIS structure functions and reduced cross-section
In this section we present a quantitative prediction of proton structure functions F 2 (x, Q 2 ) and longitudinal structure function F L (x, Q 2 ) as an outcome of our solution to the KC improved MD-BFKL equation. At small-x, the sea quark distribution is driven by gluons via. g → qq. This component from sea quark distribution can be calculated in perturbative QCD. The relevant diagram for this QCD process is shown in Fig. 11 and the contribution to the transverse and longitudinal components of the structure functions can be written using the k T -factorization theorem [57,58] where i= T, L and x x ′ is the fractional momentum carried by gluon which splits into qq pair. F (0) i includes both the quark box and crossbox contribution which comes from virtual gluon-virtual photon subprocess leading to qq production (Fig. 11 ). The gluon distribution f ( x x ′ , k 2 T ) in (4.1) represents the sum of the gluon ladder diagrams in the lower part of the box as shown in Fig. 11 is given by BFKL equation [35]. Here we will study the effect of gluon shadowing by using the function f (x, k 2 T ) which is the solution of KC improved MD-BFKL.
The explicit expression for quark box contribution F 0 i (x ′ , k 2 T , Q 2 ) can be obtained by writing four momentum in terms of the convenient light like momenta p and q ′ ≡ q + xp, where x = Q 2 2p.q and Q 2 = −q 2 like as usual (see Fig. 11 ). Now we can decompose quark and gluon momentum in terms of sudakov variables The integration should be performed over the box diagram subject to quark mass-shell constraints [38] ( which leads to the box contribution [59] F (0) where the denominators L i are Note thatF i.e. the x ′ integration of (4.1) is implicit in the d 2 k ′ T and dβ integration in (4.4) and (4.5). Now plugging (4.4) and (4.5) in the k T factorization formula (4.1) we get Equations (4.6) and (4.7) serve as the basic tool for calculating structure functions at smallx provided that the gluon distribution f (x, k 2 T ) is known. An analytical approach towards the calculation of F 2 and F L can be found in literature [60,61] for fixed coupling case using linear BFKL equation. In the literature [37], structure functions are calculated taking gluon distribution f (x, k 2 T ) from the numerical solution of unitarized BFKL equation for running coupling consideration. It has been seen that the analytical approach [61] of calculating structure functions considerably overestimates the actual (numerical) solution [37]. This is because of the fact that analytical approach neglects terms down only by powers of ln(1/z) −1 as well as it does not accommodate running of strong coupling. However, our approach is a semi-analytical one, in the sense that we take gluon distribution from our analytical solution while further integrations are performed numerically. This indeed allows us to check the feasibility of our analytical solution in describing DIS data.
In (4.1) m q denotes the quark mass and it is taken to be m q = 1.7 GeV for charm quark while massless (m q = 0) for light quarks (u, d and s). Our phenomenology is limited for light quark therefore putting m = 0 in (4.4) and replacing the quark transverse momentum .

(4.8)
After integrating over κ ′ 2 T one can arrive at T . From (4.9) it is clear thatF T . Therefore F (x ′ , k 2 T , Q 2 ) or more correctly speaking F (x ′ , k 2 T , Q 2 )/k 2 T may be considered as the structure function of an off mass shell gluon of approximate virtuality k 2 T . In [60] differential structure functions have been studied for fixed coupling and it is found that the ratio between longitudinal F L and transverse structure functionF T is 2:9 for fixed coupling approximation. We have used this ratio directly in our calculation of longitudinal structure function. Finally we have taken an assumption The advantage of taking this assumption is that we do not have to impose the possible constraint [38] coming from x/x ′ < 1 on the region of integration. Note that in principle the factorization formula (4.1) require to be run down to k 2 T = 0. The integral itself is infrared finite as both the functions f ( x x ′ , k 2 T ) and F i (x ′ , k 2 T , Q 2 ) vanish at k 2 T = 0. However, BFKL dynamics is based on perturbative QCD which is not expected to hold the nonperturbative small k 2 T physics. On the other hand, for small k 2 T the gluon distribution  Figure 12: The pQCD prediction for proton structure function F 2 obtained for two choices of shadowing: conventional (R = 5 GeV −1 ) and "hotspot" (R = 2 GeV −1 ). Data are taken from HERA (H1 [62] and ZEUS [63]) as well as fixed target experiment (NMC [64] and BCDMS [65]). Data-sets from global parametrization groups NNPDF3.1sx [27] and PDF4LHC 15 [66] is also included. The background contribution is given by vanishes linearly with the decrease in k 2 T on account of gauge invariance [37] making the contribution small. Therefore we have neglected this small contribution from small k 2 T region in our calculations of unintegrated gluon distribution f (x, k 2 T ). Before proceeding to the realistic estimation of structure functions, we add on the "background" some non-BFKL contribution, F BG i to F i since the above prediction of structure function is not enough for describing DIS data [37]. This is because (4.6) and (4.7) represent only the LLx gluon contribution and they are not the only contributions to the DIS structure functions. Although gluonic contribution is dominant at small-x, towards higher x their effect becomes weak and we cannot neglect other non-BFKL contribution. For instance, we assume that F BG i evolves like x −0.08 which is motivated from soft pomeron intercept α P (0) = 1.08 [67]. To be precise we use (4.10) An intelligent way of calculating this non-BFKL contribution is to choose x 0 at some high x and then take F BG i (x 0 , Q 2 ) from DIS data [64] which is also listed in the figure caption. The Global data-fits PDF4LHC15 Q 2 = 45 GeV 2 Figure 13: The QCD prediction for proton structure function FL obtained for two choices of shadowing: conventional (R = 5 GeV −1 ) and "hotspot" (R = 2 GeV −1 ). Data are taken from HERA H1 [68] and ZEUS [63]. Data-sets from global parametrization group PDF4LHC 15 [66] is also included. The background contribution is given by (4.10) with F BG L (x 0 = 0.1) ≈ 0.057, 0.063, 0.071 and 0.073 corresponding to four Q 2 values viz. 5 GeV 2 , 15 GeV 2 , 35 GeV 2 and 45 GeV 2 . recent HERA DIS data taken for comparison with our result can be found in [62,68] from H1 collaboration and in [63] from ZEUS collaboration. In the text [62] from H1 collaboration inclusive neutral current e ± p scattering cross-section data collected during the years (2003-2007) is presented. The beam energies E p of corresponding H1 experiment run are 920, 575 and 460 GeV 2 . Corresponding F L data of H1 is taken from [68] where measurement are perfomed at c.m.s. energy √ s = 225 and 252 GeV 2 . On the other hand, in [63] from ZEUS collaboration reduced cross-sections for ep scattering for different cms energies viz. 318, 251 and 225 GeV 2 is presented. The fixed target data from New Muon Collaboration (NMC) [64] and BCDMS [65] collaboration which exist for x > 10 −2 is also shown in the Figure 12. Finally we have included data from global DIS data analysis paramaterization groups viz. NNPDF3.1sx [27] and PDF4LHC15 [66] for comparison. Both of the LHAPDF datasets NNPDF3.1sx [27] and PDF4LHC15 [66] include HERA as well as recent LHC data with high precision PDF sensitive measurements.
The x dependence of the structure functions F L and F 2 is shown in Fig. 12 and Fig. 13 for the four Q 2 values 5 GeV 2 , 15 GeV 2 , 35 GeV 2 and 45 GeV 2 . QCD prediction shows a satisfactory agreement with DIS data for both structure functions F 2 and F L . Interestingly both data and theory for F 2 structure function seem to preserve the x −λ singular behavior rather than showing taming due to net shadowing effect in the asymptotic limit x → 0. It is difficult to observe the existence of any sizable shadowing effect for x > 10 −3 . However for x < 10 −3 , recalling that the shadowing term is proportional to 1/R 2 it is seen that for R = 2 GeV −1 (gluons in hotspots) the rise of structure function becomes slower than that of R = 5 GeV −1 (R = R H , hadron radius) as expected. But even the extreme form of the gluon shadowing (R = 2 GeV −1 ) suppresses F 2 only by about 10% or less at 10 −3 . It just depicts that shadowing has the negligible impact on structure functions even towards smaller x. This is because of the fact that in this low x regime gluons are expected to drive the sea quark distribution via g → qq. Therefore a similar sea quark contribution to the structure function F 2 in addition to the gluon contribution can be expected in low x regime. On the other hand from Fig. 13 it is clear that for x > 10 −3 , the size of the longitudinal structure function is negligibly small while for very small-x regime (x < 10 −3 ) , F L grows eventually. This is in accord with our expectation since measurement of F L directly probes the gluonic content of the proton which is dominant in small-x regime. The proton structure functions F 2 and F L are of complementary in nature. These are related to the γ * p cross-sections of longitudinally and transversely polarized photons σ L and σ T as F L ∝ σ L and F 2 ∝ (σ L + σ T ). Since σ L and σ T are positive, that imposes a restriction on F L and F 2 i.e. 0 ≤ F L ≤ F 2 . To circumvent the need of a proper relationship between the structure functions and polarized cross-sections to study γ * p cross-section one can define a ratio between the structure functions or the the equivalent cross-section ratio (4.11) The advantage behind this formulation is that this is independent of any normalization factors.
In Fig 14 a phenomenological comparison between data and theory is shown to illustrate Q 2 dependence of the ratio R(x, Q 2 ). The HERA data taken for comparison can be found in [69] where ep cross-section data measured for two center of mass energies √ s = 225 and 252 GeV is enlisted. The corresponding center of mass energy √ s and inelasticity y for measurement of R is also given in the figure caption. From Fig. 14 it is clear that the HERA data is in general agreement with the QCD expectation. Available H1 data for very low Q 2 (Q 2 ≤ 5 GeV) is also sketched in Fig. 14. We have excluded the very low Q 2 region in our calculation of R since we are bounded to stick in the perturbative region. However, a more realistic study in this transition region is performed in Sect. 4.2. Now we try to predict reduced cross-section σ r for ep scattering process from the knowledge of structure functions that we have discussed above. The inclusive deep-inelastic differential cross-section for ep scattering can be represented in terms of three structure functions F 2 , F L and xF 3 . The structure functions have direct dependence with DIS kinematic variables, x and Q 2 only, while the cross-section has additional dependence with inelasticity y = Q 2 /sx. The inclusive cross-section for neutral current ep scattering is given by where Y ± = 1±(1−y) 2 . The cross-section for exchanged virtual boson (Z or γ * ) is related to F 2 and F 3 in which contribution from both longitudinal and transverse boson polarization state exists. On the other hand, only the longitudinal polarized virtual boson exchange processes contribute to F L which has a significant impact on higher order QCD though it vanishes at lowest order. At small momentum transfer Q 2 (i.e. Q 2 ≪ M 2 Z ≈ 800 GeV 2 ), interaction of massless photon is dominant over the exchange of heavy Z boson. Thus the contribution from Z boson exchange and the influence of the term xF 3 is negligible at low and moderate Q 2 . Therefore in this range of Q 2 , in one photon exchange approximation, the differential cross-section formula (4.12) can be written as where σ r is the reduced cross-section. Note that (4.13) is symmetric under charge exchange i.e. identical for both e + p and e − p processes unlike (4.12). Additionally, it is independent of incoming electron helicity state. Thus at low and moderate Q 2 (≤ m 2 Z ) which is our region of study, the knowledge of F 2 and F L is enough to predict the reduced cross-section. We can also express reduced cross-section in terms of the ratio R(x, Q 2 ) defined in (4.11) replacing F L by F L = R 1+R F 2 which yields (4.14) The x dependence of e ± p reduced cross-section σ r calculated for center of mass energy √ s = 318 GeV is shown in Fig. 15. Our theoretical expectation is compared with HERA H1 measurement [70,71] and ZEUS ( √ s = 318 GeV) [63]. The available H1 data for low Q 2 (≤ 12 GeV 2 ) measured at √ s = 318 GeV (SVX) and 300 GeV (NVX-BST) is taken from [70], while the same for high Q 2 (> 12 GeV 2 ) is taken from [71]. The cross-section measurement of SVX is found to be slightly higher than that of NVX-BST as expected because of the increase in center of mass energy. The error band included in data represents the total uncertainty determined by the quadratic sum of systematic and statistical uncertainties. Both theory and data agree well with each other for our phenomenology range Q 2 ≤ 100 GeV 2 . For each Q 2 , starting at some high x the reduced cross-section first increases as x → 0 and then an abrupt fall in cross-section can be observed in both theory and data at very small-x region (x < 10 −4 ) . For all Q 2 , this region of x corresponds to the highest inelasticity y ≈ 0.65 (y = Q 2 /sx) and thus suppression of cross-section at y ≈ 0.65 can be attributed to the influence of F L . In simple words, towards very small-x (or high y) the monotonic rise of F 2 is suppressed by the contribution from longitudinal structure function F L thereby causing an overall fall in the cross-section. For low inelasticity y < 0.65, the contribution from the longitudinal structure function is small on the other hand structure function F 2 exhibits a steady increase as x → 0. Therefore in the region where x is not so small, the growth of the cross-section is found to be power like as expected.

Virtual photon-proton cross-section prediction in transition region
Generally the photon-proton interaction is classified into two separate processes depending upon the photon virtuality Q 2 viz. photoproduction (at low Q 2 ) and deep inelastic scattering (at high Q 2 ). DIS is considered as a basic tool in pQCD where the point like virtual photon directly probes the partonic contents of the proton. On the other hand, photoproduction is completely nonperturbative phenomena defined in the limit of vanishing Q 2 where real (or quasi-real) photons interact with the proton more likely a hadron-hadron collision. The experiment at HERA collider provides a unique opportunity to study both the processes photoproduction and DIS on their respective kinematic domains. The Q 2 dependence of the proton structure functions are well described by pQCD over a wide range of x and Q 2 [72,73] in accordance with HERA data. However, for Q 2 2 GeV 2 (photo production region) the pQCD breakdowns since the higher order contributions to the perturbative expansion becomes very large. In this region, data can be only described by non perturbative phenomenological models [74]. Our present study is especially focused on the transition region (2 GeV 2 ≤ Q 2 ≤ 10 GeV 2 ) from photoproduction to deep inelastic scattering. To measure the photon-proton cross-section in the transition region two dedicated runs were performed in the years 1999 (Nominal vertex "NVX'99") and 2000 (Shifted Vertex "SVX'00") by H1 experiment at HERA. The published data can be found in [70].
Recall the neutral current ep double differential cross-section formula (4.13) defined in the region Q 2 ≪ M 2 Z . In this region massive boson (M Z ) exchange is neglected and only one photon exchange is considered, thereby the role of incoming electron reduces to be a source of virtual photon interacting proton. Thus we can recast the formula (4.13) in terms of photon-proton reaction. In fact the structure function F 2 and F L in (4.13) related to the longitudinally and transversely polarized polarized photon-proton scattering cross-sections σ L and σ T by the relations which hold good at low x. Considering (4.15) and (4.16) the reduce cross-section in (4.13) can be written as where is the effective virtual photon-proton cross-section. Note that expression for σ eff γ * p is similar to the total crosssection σ tot γ * p which is linear combination of σ L and σ T i.e. σ tot γ * p = σ L + σ T . In fact the total cross-section σ tot γ * p and σ eff γ * p can be regarded as the same quantity at low inelasticity y i.e. σ eff γ * p y 2 →0 −−−→ σ tot γ * p which differ only in the region of high y. Figure 16 shows the measurement of the virtual photon-proton cross-section σ * γ * p as a function of Q 2 corresponding to different values of W . The total cross-section σ tot γ * p is often expressed as a function of Q 2 and invariant mass W . The standard relation between W , x and Q 2 is W = Q 2 (1 − x)/x. Since Q 2 ≈ syx, for small-x we can have an approximate relationship between W and y i.e. W 2 ≃ sy which we have used in our calculations. HERA included for comparison with our theoretical prediction. We have chosen √ s = 318 GeV for our calculation analogous to H1 measurement. Both the HERA data and theoretically measured cross-sections for different values of W are multiplied by the factors multiple of 2 as indicated in the figure. A good agreement between the theory and HERA data can be observed in Fig. 16 for the wide range of W . This analysis ensures the compatibility of our theory in low Q 2 transition region from photoproduction to deep inelastic scattering.

Extraction of colinear gluon density from F L
In this section, we study the role of longitudinal structure function F L in the measurement of colinear gluon distribution, xg(x, Q 2 ). Approximate relationship between F L and xg(x, Q 2 ) can be found in literature [76][77][78] which is often used in extraction of xg(x, Q 2 ) from F L by various collaboration and global DIS data parameterization groups.
Before going to the calculation part we want to revisit the basic physics behind the correlation between longitudinal structure function and gluon distribution. In QCD the polarised structure functions F L and F T can be represented as the linear combination of F 1 and F 2 structure functions: In naive parton model which is based only on charged massless spin 1 2 partons with no gluon, the longitudinal structure function F L vanishes. That leads to nothing but the x . 10 4 Figure 17: Q 2 dependence of colinear gluon density xg(x, Q 2 ) extracted from longitudinal structure function F L using (4.21). The theoretical prediction from KC improved MD-BFKL is compared with HERA H1 [68] data.
Callan-Gross relation F 2 − 2xF 1 = 0. But in QCD where the gluons come into the picture, higher order QCD corrections change the Calan-Gross relation leading to (4.20) and F L becomes proportional to α s . Thus measurement of F L directly probes the gluonic content of the proton at small-x. However, it should be noted that the contribution from F L to the cross-section is sizeable only when y is large (see (4.13)) i.e. either at small-x or at large Q 2 (since y ∝ Q 2 x ). In [79] Altarelli-Martinellin have proposed a direct relationship between gluon distribution xg(x, Q 2 ) and F L . But the difficulty in this relation is that it can't be solved analytically. However an approximate solution can be found in [76] xg(x, where we have neglected the tiny contribution from F 2 . Equation (4.21) suggests that gluon density can be simply read off measurement of F L by rescaling the value of x by 0.4. We use this relation to measure colinear gluon density from the knowledge of longitudinal structure function studied in Sect. 4.1.
The Q 2 dependence of Gluon density xg(x, Q 2 ) averaged over x is shown in Fig. 17. Our result is compared with that of H1 collaboration [68] in the kinematic range 1.5 GeV 2 ≤ Q 2 ≤ 800 GeV 2 . A nonlinear fit to both data and theory is also shown. Our theoretical prediction shows a rough agreement with HERA data. Towards large Q 2 the growth of our gluon distribution prediction looks less satisfactory compared to the DIS data. This may be because of α s (Q 2 ) in (4.21), for which the variation of xg(x, Q 2 ) with respect to Q 2 is not rapid. Moreover, our theory is only compatible for low and moderate Q 2 (Q 2 ≪ M 2 Z ) in the  Figure 18: Small-x dependence of colinear gluon density extracted from longitudinal structure function F L using (4.21). The theoretical prediction from KC improved MD-BFKL is compared with that of NNLO solution of GLR-MQ equation [80,81] as well as various global datasets viz. NNPDF3.1sx [27], HERAPDF 15 [28], ABMP [29], ATLAS [30] and CT14 [31].
approximation of single photon exchange and it is not an appropriate choice for study DIS data in large Q 2 domain where production of heavy gauge boson M Z can't be neglected.
In Fig. 18 the x evolution of colinear gluon density extracted according to (4.21) is compared to the various global data analysis parameterization groups viz. NNPDF3.1sx [27], HERAPDF 15 [28], ABMP [29], ATLAS [30] and CT14 [31]. Additionally, we have included prediction from nonlinear GLR-MQ equation which is solved analytically up to NNLO accuracy in [80,81]. We have plotted the gluon density in the kinematic range of 10 −6 ≤ x ≤ 10 −2 at four fixed values of photon virtuality Q 2 viz., 10, 35, 50 and 100 GeV 2 respectively. From the graph, it can be seen that our predicted gluon density rises towards decreasing x, which is in agreement with perturbative QCD fits at small-x. The gluon density obtained from (4.21) is seemed to rise slowly as compared to the predicted gluon density of the other PDF groups as shown in the figure, this may be attributed to the effect of nonlinearity in the equation. However, our theoretical result predicts slightly higher gluon density than that of the other PDF groups in the region of x ≥ 10 −3 at Q 2 = 50 and 100 GeV 2 , as is visible in the graph. can be found in the literature [45]. Surprisingly this parametrization provides a very good description of F 2 DIS data for unified BFKL-DGLAP framework. Now recalling and assuming that f (x, k 2 T ) ∼ k 4 for low k 2 T (i.e. k 2 T < k 2 0 ) [85] we obtain Our result is compared with that of GBW parametrization. We set the r evolution range as 0 ≤ r ≤ 2 GeV −1 considering small dipole approximation compared to transverse size of the target. It is observed that our prediction of dipole crosssection shows a similar behavior as GBW parametrization. However, the small disagreement seen is probably because of the difference in the choices of input distribution in the two approaches while calculating σ(r, x).

Conclusion
In conclusion, we have presented a phenomenological study on the behavior of unintegrated gluon distribution at small-x and moderate k 2 T region particularly 10 −5 ≤ x ≤ 10 −2 and 2 GeV 2 ≤ k 2 T ≤ 100 GeV 2 which is also the accessible kinematic range to experiments performed at HERA ep collider. In the beginning, we have improved the MD-BFKL equation supplementing so-called kinematic constraint on it. Then we have solved this unitarized BFKL equation analytically in order to study x and k 2 T dependence of unintegrated gluon distribution function. We have also performed an indirect comparison of our prediction with TMDlib datasets viz. PB-NLO-HERAI+II-2018 [24,25] and ccfm-set-A0 [26]. The x evolution of gluon distribution shows the singular x −λ type behavior of gluon evolution tamed by shadowing correction. We found that for intense shadowing (R = 2 GeV −1 ), towards smaller x, certainly from x = 10 −3 , the gluon distribution emerges from rapid BFKL growth which indicates the dominance of gluon shadowing. While in case of conventional shadowing (R = 5 GeV −1 ) an appreciable modification of BFKL behavior can only be seen from x = 10 −4 . The k 2 T dependence of unintegrated gluon distribution is also studied. Although obvious suppression due to net shadowing correction is seen in our studies however no significant saturation phenomenon is observed. The reason is the nonlinear contribution term is suppressed by the factor 1/k 2 T at large values of k 2 T . We have obtained a more general solution of KC improved MD-BFKL equation implementing a pre-defined input distribution on it. This has allowed us to visualize the gluon evolution in three dimensions R 3 into x-k 2 T phase space. We have also shown the sensitiveness of unintegrated gluon distribution towards R and λ using density plots in x-k 2 T plane.
A remarkable achievement of obtaining an analytical solution in this work is its implication on qualitative studies on geometrical scaling which is presented in Sect. 3. Starting from a basic concept of differential geometry and knowledge of level curves we have managed to obtain an equation of the critical boundary which is supposed to separate low and high gluon density regions.
In Sect. 4 we have studied the small-x dependence of the structure function F 2 and F L obtained via k T -factorization formula F i = f ⊗ F (0) i . Here the unintegrated gluon distribution f (x, k 2 T ) is taken from our analytical solution of KC improved MD-BFKL equation setting boundary condition at x 0 = 0.01. Surprisingly the quantitative size of the shadowing correction to F 2 and F L is found to be very small and the structure functions seem to hold the singular x −λ behavior. Even at intense shadowing condition (R = 2 GeV −1 ) the F 2 structure function is found to be suppressed only by 10%. In addition we have measured e ± p reduced cross-section as well as equivalent γ * p longitudinal to transverse cross-section ratio R(x, Q 2 ) from the knowledge of F 2 and F L . Our results are compared with recent high precision HERA measurements. The comparison reveals a good agreement between our theory and DIS data.
In Sect. 4.2 we have examined the virtual photon-proton effective cross-section, particularly in the transition region from photoproduction to deep inelastic scattering. The quantity σ eff γ * p serves the role of the total cross-section if small inelasticity y is concerned. We have compared our predicted σ eff γ * p with the HERA data of two dedicated runs "NVX'99" and "SVX'00" by H1 as well as ZEUS for the transition region. Our theoretical prediction shows well-consistency with HERA data in the transition region.
In Sect. 4.3 we have presented a phenomenological measurement of colinear gluon density xg(x, Q 2 ) from longitudinal structure function, F L using the approximate relation (4.21). For comparison, we have plotted NNLO solution of GLR-MQ equation as well as various global data fits viz. NNPDF3.1sx, HERAPDF 15, ABMP, ATLAS, and CT14. The Q 2 dependence of the colinear gluon density is also investigated and a satisfactory agreement with HERA H1 data is observed. Finally in Sect. 5 we have presented a dipole picture of our KC improved MD-BFKL equation. In this dipole phenomenology, we have shown that a simple transformation of original MD-BFKL equation leads to the BK equation whereas the same prescription leads KC improved MD-BFKL equation to KC improved BK equation near the saturation limit. In this section, we have also calculated dipole cross-section σ(r, x) and sketch a comparison with that of GBW parametrization.
In summary, there are several attractive features of our present study. First, we have able to predict a wide range of physical quantities, starting right from our solution for unintegrated gluon density. Secondly, all analysis is performed in terms of relatively small numbers of parameters. Moreover, the idea developed in this work for studying geometrical scaling can be implemented in any other framework. Finally, a very significant feature of this analysis is that we have considered two extreme possibilities of shadowing which can be distinguished by experimental data. In the end, we conclude that as per feasibility towards HERA DIS data is concerned, the KC improved MD-BFKL equation can be considered as a reliable framework for exploring high energy physics over a wide range of x and k 2 T which is also relevant for LHC probe and future collider phenomenology.