∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>
11institutetext: Department of Physics, Tezpur University, Sonitpur, Assam,
India

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

P. Phukane1,addr1  M. Lalunge2,addr1  J. K. Sarmae3,addr1

###### Abstract

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.

††journal: Eur. Phys. J. C

###### Contents

  1. 1 Introduction
  2. 2 Nonlinear effects in kinematic constraint improved MD-BFKL equation
    1. 2.1 Construction of kinematic constraint improved MD-BFKL
    2. 2.2 Analytical solution of KC improved MD-BFKL
      1. 2.2.1 $x$ and $k_{T}^{2}$ evolution
      2. 2.2.2 Complete solution of KC improved MD-BFKL
  3. 3 Equation of Critical line and prediction of saturation scale: a differential geometric approach
  4. 4 KC improved MD-BFKL prediction of HERA DIS data
    1. 4.1 DIS structure functions and reduced cross section
    2. 4.2 Virtual photon-proton cross section prediction in transition region
  5. 5 Conclusion

##  1 Introduction

In perturbative QCD the two well-known parton evolution equations DGLAP 63 ;
64 ; 65  and BFKL 66 ; 67  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 $\alpha_{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 $(\alpha_{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 68 . 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\rightarrow g$) which is simply the
inverse process of gluon splitting ($g\rightarrow g+g$). The modification of
gluon recombination to the original DGLAP equation was first performed in GLR-
MQ equation 16 ; 75 ; 86 ; 87  proposed by Gribov-Levin-Ryskin and Mueller-
Qui. 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 59 ; 69 . The GLR-MQ equation was originally derived
using AGK-cutting rules 71  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 70 . 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 69  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 collinear
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 1  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 10 ; 11  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 85 . 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, like GLR-MQ equation the BK equation is also irrelevant to
the gluon antishadowing. On the other hand, besides inclusion of antishadowing
correction term, the nonlinear terms in MD-BFKL equation hold a 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 1 . In the literature 73  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:

$\displaystyle q_{T}^{2}$ $\displaystyle<\frac{k_{T}^{2}}{z}\text{ }\text{ LDC
\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{6,74}{\@@citephrase{(}}{\@@citephrase{)}}}},$ (1)
$\displaystyle q_{T}^{2}$ $\displaystyle<\frac{(1-z)k_{T}^{2}}{z}\text{
\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{3}{\@@citephrase{(}}{\@@citephrase{)}}}},$ (2)
$\displaystyle k_{T}^{{}^{\prime}2}$ $\displaystyle<\frac{k_{T}^{2}}{z}\text{
}\text{ BFKL \cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{3}{\@@citephrase{(}}{\@@citephrase{)}}}}.$ (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^{{}^{\prime}2}\approx k_{T}^{{}^{\prime}2}$. The kinematic
constraint gives an implicit cutoff on $k_{T}^{{}^{\prime}2}$ as depicted by
(3). The inequality (3) can be considered as a special case of (1) recalling
the fact that for a given value of $k_{T}^{2}$, a high $q_{T}^{2}$ implies an
equally high $k_{T}^{{}^{\prime}2}$ 3 . Although there are other constraints
coming from energy-momentum conservation 7 , this bound is considerably
tighter than the later 3 . 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_{T}^{2}$ evolution followed by three-
dimensional realization of gluon evolution in $x\text{-}k_{T}^{2}$ phase space
(Sect. 2.2). We have also extracted collinear gluon distribution $xg(x,Q^{2})$
from unintegrated gluon distribution $f(x,k_{T}^{2})$ and compared our
prediction with that of global parameterization groups NNPDF 3.1sx and CT 14.
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. Finally, in Sect. 5 we summarize and outline our conclusion as well
as possible future prospects.

##  2 Nonlinear effects in kinematic constraint improved MD-BFKL equation

###  2.1 Construction of kinematic constraint improved MD-BFKL

The modified BFKL (MD-BFKL) equation reads 1 ,

$\begin{split}&-x\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial x}=\\\
&\frac{\alpha_{s}N_{c}k_{T}^{2}}{\pi}\int_{k_{T_{\min}}^{{}^{\prime}2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left[\frac{f(x,k_{T}^{{}^{\prime}2})-f\left(x,k_{T}^{2}\right)}{|k_{T}^{{}^{\prime}2}-k_{T}^{2}|}+\frac{f\left(x,k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\right]\\\
&-\frac{36\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}(x,k_{T}^{2})+\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}(\frac{x}{2},k_{T}^{2}),\end{split}$
(4)

where $f\left(x,k_{T}^{2}\right)$ denotes the gluon distribution unintegrated
over the transverse momentum of gluon $k_{T}$ and
$k_{T_{\min}}^{{}^{\prime}2}$ is the infrared cutoff of the evolution. The
first part of (4) linear in $f(x,k_{T}^{{}^{\prime}2})$ (or $f(x,k_{T}^{2})$)
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_{T}^{{}^{\prime}2})$ 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_{T}^{{}^{\prime}2}=k_{T}^{2}$ cancels between real and virtual
contributions.

The quadratic terms in (4) viz. $\frac{\partial
f_{\text{shad}}}{\partial\ln{1/x}}=-\frac{36\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}$ and $\frac{\partial
f_{\text{antishad}}}{\partial\ln{1/x}}$ $=\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}$ depict shadowing and
antishadowing correction to the original BFKL equation respectively. The
factor $\pi 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\approx
R_{H}=5\text{ 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\text{ GeV}^{-1}$. In the latter case $(R=2\text{
GeV}^{-1})$ we have ”hot spots” 75 ; 72 .

Figure 1:  From left: (a) chain of sequential gluon emission which forms the
basis of BFKL equation. On squaring the amplitude of (a) the ladder diagram
(b) is generated which when summed gives the BFKL kernel. (c) and (d) are two
simple examples of the inhomogeneous driving term of (10) which correspond to
the shaded region of gluon-virtual photon coupling in (b).

The interpretation of the real emission term in the BFKL kernel is that we
have a splitting $k^{{}^{\prime}}\rightarrow k+q$ inside hadron resulting an
infinite chain of reggeized gluons labeled as
$k_{1},k_{2},k_{3},\text{...}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.

$\displaystyle x_{1}\ll x_{2}\ll x_{3}\text{.....}\ll x_{n-1}\ll x_{n},$ (5)
$\displaystyle k_{1T}\sim k_{2T}\sim k_{3T}\sim\text{...}.\sim k_{n-1T}\sim
k_{\text{nT}}.$ (6)

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,

$\displaystyle k^{2}=2k^{+}k^{-}-k_{T}^{2}\approx k_{T}^{2}.$ (7)

The above kinematics corresponding to (5),(6) and (7) is referred as multi-
Regge kinematics. The inequality (5) implies
$z=\frac{x}{\frac{x}{z}}=\frac{x_{n-1}}{x_{n}}\ll 1$ and since transverse
momenta are of the same order: $k_{T}\simeq k_{T}{}^{{}^{\prime}}$, depict a
cutoff

$k_{T}{}^{{}^{\prime}}{}^{2}<\frac{k_{T}{}^{2}}{z},$ (8)

which is so-called kinematic constraint or consistency constraint for real
gluon emission 3 ; 4 ; 5 ; 6 . The BFKL kernel in (4) is at leading
logarithmic $1/x$ (LLx) accuracy. However, higher order corrections to the
BFKL equation are already been evaluated up to NLLx accuracy, which turns out
to be quite large 2 . Implementation of the constraint (8) 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 NLLx correction in BFKL
kernel.

Recalling that BFKL equation can be written as an integral equation of
$f(x,k_{T}^{2})$ 15 ; 16 , the kinematic constraint (8) can be imposed onto
the real emission part of the BFKL kernel as follows

$\begin{split}\int_{x}^{1}\frac{dz}{z}\int_{k_{T_{\min}}^{{}^{\prime}2}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\bigg{[}&\frac{\Theta\left(\frac{k_{T}^{2}}{z}-k_{T}^{{}^{\prime}2}\right)f\left(\frac{x}{z},k_{T}^{{}^{\prime}2}\right)-f\left(\frac{x}{z},k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}}{}^{2}-k_{T}^{2}\right|}\\\
&+\frac{f\left(\frac{x}{z},k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]},\end{split}$
(9)

where the heaviside function
$\Theta\left(\frac{k^{2}}{z}-k^{{}^{\prime}}{}^{2}\right)$ in (9) serves the
purpose of the kinematic cutoff. The upper limit of the integration over
$k_{T}^{{}^{\prime}2}$ is implicit in $\Theta$. The BFKL kernel gains one more
degree of freedom after implementation of kinematic constraint for real
emission. Now expressing the MD-BFKL equation (4) in terms of the KC improved
BFKL kernel (9) yields

$\begin{split}&f\left(x,k_{T}^{2}\right)=f^{(0)}\left(x,k_{T}^{2}\right)+\\\
&\frac{\alpha_{s}k_{T}^{2}N_{c}}{\pi}{\int_{x}^{1}\frac{dz}{z}\int_{k_{T_{\min}}^{{}^{\prime}2}}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\bigg{[}\frac{\Theta\left(\frac{k_{T}^{2}}{z}-k_{T}^{{}^{\prime}2}\right)f\left(\frac{x}{z},k_{T}^{{}^{\prime}2}\right)-f\left(\frac{x}{z},k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&+\frac{f\left(\frac{x}{z},k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]}-\frac{36\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}\int_{x}^{1}\frac{dz}{z}f^{2}\left(x,k_{T}^{2}\right)\\\
&+\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}\int_{x}^{1}\frac{dz}{z}f^{2}\left(\frac{x}{2},k_{T}^{2}\right).\end{split}$
(10)

The above equation (10) is the integral form of our KC improved MD-BFKL
equation. The inhomogeneous driving term $f^{(0)}(x,k_{T}^{2})$ depicts gluon-
proton coupling corresponding to the shaded region in Fig. 1(b). The Fig. 1(c,
d) represents two simple possible contributions to $f^{0}(x,k_{T}^{2})$ which
indicates radiation of gluons from valence quark.

To sketch an integro-differential form of KC improved MD-BFKL equation out of
(10) we differentiate (10) w.r.t. $\ln(1/x)$, then using properties of
$\Theta$ and Dirac-$\delta$ function, to be specific
$\Theta^{\prime}(t)=\delta(t)$ and $f(t)\delta(t-a)=f(a)\delta(t-a)$ and doing
some algebra one can show

$\begin{split}\frac{\partial}{\partial\ln\frac{1}{x}}\int_{x}^{1}\frac{\text{d}z}{z}\Theta\left(\frac{k_{T}^{2}}{z}-k_{T}^{{}^{\prime}2}\right)&f(x,k_{T}^{{}^{\prime}2})\longrightarrow\Theta(k_{T}^{2}-k_{T}^{{}^{\prime}2})f(x,k_{T}^{{}^{\prime}2})\\\
&+\Theta(k_{T}^{{}^{\prime}2}-k_{T}^{2})f\left(\frac{k_{T}^{{}^{\prime}2}}{k_{T}^{2}}x,k_{T}^{{}^{\prime}2}\right).\end{split}$
(11)

Above prescription allows one to express (10) in the following integro-
differential form:

$\displaystyle\begin{split}&-x\frac{\partial
f\left(x,k_{T}^{2}\right)}{\partial
x}=\frac{\alpha_{s}k_{T}^{2}N_{c}}{\pi}{\int_{k_{T_{\min}}^{{}^{\prime}2}}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\\\
&\times\bigg{[}\frac{\Theta\left(k_{T}^{2}-k_{T}^{{}^{\prime}2}\right)f\left(x,k_{T}^{{}^{\prime}2}\right)+\Theta(k_{T}^{{}^{\prime}2}-k_{T}^{2})f\left(\frac{k_{T}^{{}^{\prime}2}}{k_{T}^{2}}x,k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&-\frac{f(x,k_{T}^{2})}{|k_{T}^{{}^{\prime}2}-k_{T}^{2}|}+\frac{f\left(x,k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]}-\frac{36\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}\left(x,k_{T}^{2}\right)\\\
&+\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}\left(\frac{x}{2},k_{T}^{2}\right).\end{split}$
(12)

In derivation of (12) we have neglected the term $x\frac{\partial
f^{(0)}}{\partial x}$. This is justified in the sense that $x\frac{\partial
f^{(0)}}{\partial x}$ is much less singular than $x\frac{\partial f}{\partial
x}$ at small-$x$. Moreover if we consider $f^{(0)}$ is particularly
contributed by the diagrams shown in Fig. 1 (c, d) then $f^{(0)}$ would be
independent of $x$ (or, $\frac{\partial f^{(0)}}{\partial x}=0$ ) at small-$x$
limit 6 ; 15 .

Figure 2:  QCD cut diagram representing nonlinear terms in (4) with
$4\rightarrow 2$ gluon recombination kernel (dashed box). The dashed box
includes all possible pQCD diagrams which couple 4 gluons to 2 gluons.

To simplify the distribution function
$f(\frac{k_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}x,k_{T}^{2})$
corresponding to real emission term in (12) 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 (${\alpha_{P}(0)\approx
1.08}$) and Reggons with intercept $\alpha_{R}(0)\approx 0.5$. The Regge
behavior corresponding to sea quark and anti-quark distributions is given by
$q_{sea}\sim x^{-\alpha_{P}}$ whereas that of valence quark distribution is
given by $q_{v}\sim x^{-\alpha_{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

$\alpha_{P}^{BFKL}(0)=1+\lambda_{\text{BFKL}},$

where $\lambda_{\text{BFKL}}=\frac{3\alpha_{s}}{\pi}28\zeta(3)$, $\zeta$ being
Reiman zeta function 9 . The typical value of $\lambda_{\text{BFKL}}$ for
$\alpha_{s}$=0.2 is $\sim 0.5$ implies $\alpha_{P}^{BFKL}(0)\approx 1.5$ which
is potentially large in magnitude compared to soft pomeron
$(\alpha_{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^{-\lambda_{i}}$ (i=$\sum$ (singlet structure
function)and g (gluon distribution)), where$\lambda_{i}$ is the pomeron
intercept minus one $(\alpha_{P}(0)-1)$. At small-$x$, the leading order
calculations in $\text{ln}(1/x)$ with fixed $\alpha_{s}$ predicts a steep
power law behavior of $\textit{f}\left(x,k^{2}\right)\sim
x^{-\lambda_{\text{BFKL}}}$ 7 ; 8 . This motivates us to consider a simple
form of Regge factorization as follows,

$\begin{split}f\left(\frac{k_{T}^{{}^{\prime}2}}{k_{T}^{2}}x,k_{T}^{2}\right)&\simeq
x^{-\lambda_{\text{BFKL}}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda_{\text{BFKL}}}H(k_{T}^{2})\\\
&=\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}f\left(x,k_{T}^{2}\right),\end{split}$
(13)

where we have dropped the subscript on $\lambda$ and in the rest of the text
we will follow this notation. Similar way we can express

$\begin{split}f\left(\frac{x}{2},k^{2}\right)&\simeq
2^{\lambda}f\left(x,k^{2}\right).\end{split}$ (14)

This type of Regge like form considered in (13) and (14) is supported by
various literatures 76 ; 77 ; 78 ; 79 . But note that Regge factorization
can’t be taken as good ansatz for the entire kinematic region of $x$ and
$k_{T}^{2}$ 80 . This type of Regge behavior is considered to be valid only in
the vicinity of the saturation scale, $Q_{s}^{2}$ where scattering amplitude
depends only on a single dimensionless variable, $Q^{2}/Q_{s}^{2}$. The Regge
theory is applicable if the quantity invariant mass, $W$
$(=\sqrt{Q^{2}(1-x)/x})$ is much greater than all other variables. Therefore,
we expect it to be valid if $x$ is enough small, for any value of $Q^{2}$.

Now recalling $\Theta(t)=1-\Theta(-t)$ and substituting (13) and (14) in (12)
we get

$\begin{split}&-x\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial x}=\\\
&\frac{\alpha_{s}N_{c}k_{T}^{2}}{\pi}\int_{k_{T_{\min}}^{{}^{\prime}2}}\frac{dk^{\prime
2}}{k_{T}^{{}^{\prime}}{}^{2}}\bigg{[}\Theta\left(k_{T}^{2}-k_{T}^{{}^{\prime}2}\right)\frac{1-\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}f\left(x,k_{T}^{{}^{\prime}2}\right)\\\
&+\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{f\left(x,k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}2}-k^{2}\right|}-\frac{f\left(x,k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}2}-k^{2}\right|}+\frac{f\left(x,k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]}\\\
&-\frac{36\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}(1-2^{2\lambda-1})f^{2}\left(x,k_{T}^{2}\right).\end{split}$
(15)

The above equation (15) 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 (15). As we follow, in the region far below
saturation limit, for the canonical choice of $\lambda\sim 0.5$, the quadratic
term in (15) tends to vanish since $(1-2^{2\lambda-1})\rightarrow 0$. This
indicates that below saturation region both the contribution from shadowing
correction $\frac{\partial f_{\text{shad}}}{\partial\ln{\frac{1}{x}}}$ and
antishadowing correction $\frac{\partial
f_{\text{antishad}}}{\partial\ln{\frac{1}{x}}}$ coexists but seem to balance
each other for which nonlinear effect becomes negligible. Therefore, in the
below saturation regime, we can revert (15) back to the original BFKL equation
choosing $\lambda=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 (14) invalid. Rather essentially, we should take the
approximation

$f\left(\frac{x}{2},k_{T}^{2}\right)\simeq f(x,k_{T}^{2}),$

which is supposed to justify our understanding of the saturation region.
Taking this approximation in (12) one can arrive at

$\begin{split}&-x\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial x}\\\
&=\frac{N_{c}\alpha_{s}k_{T}^{2}}{\pi}\int_{k_{T_{\min}}^{{}^{\prime}2}}\frac{dk^{\prime
2}}{k_{T}^{{}^{\prime}}{}^{2}}\bigg{[}\Theta\left(k_{T}^{2}-k_{T}^{{}^{\prime}2}\right)\frac{1-\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}f\left(x,k_{T}^{{}^{\prime}2}\right)\\\
&+\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{f\left(x,k_{T}^{{}^{\prime}2}\right)}{\left|k_{T}^{{}^{\prime}2}-k^{2}\right|}-\frac{f\left(x,k_{T}^{2}\right)}{\left|k_{T}^{{}^{\prime}2}-k^{2}\right|}+\frac{f\left(x,k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]}\\\
&-\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}f^{2}\left(x,k_{T}^{2}\right).\end{split}$
(16)

Equation (16) serves the evolution of gluon near the saturation limit. Here
the $\lambda$ dependence of the nonlinear term has been completely wiped up,
unlike (15). In this regime, the shadowing contribution becomes twice as that
of antishadowing effect forecasting a net shadowing effect in the evolution.

###  2.2 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
$\alpha_{s}$, so do fix $\lambda$, therefore, the solution for both (15) and
(16) 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}^{{}^{\prime}}\backsimeq k_{T}$ in BFKL multi-Regge kinematics, this
allows us to write the gluon distribution in Taylor series,

$f\left(x,k_{T}^{{}^{\prime}2}\right)=f\left(x,k_{T}^{2}\right)+\frac{\partial
f\left(x,k_{T}^{2}\right)}{\partial
k_{T}^{2}}\left(k_{T}^{{}^{\prime}2}-k_{T}^{2}\right)+\mathcal{O}\left(k_{T}^{{}^{\prime}2}-k_{T}^{2}\right),$
(17)

where $\mathcal{O}\left(k_{T}^{{}^{\prime}2}-k_{T}^{2}\right)$ denotes the
higher order terms. Above series is a convergent series in
$\left(k_{T}^{{}^{\prime}2}-k_{T}^{2}\right)$ as no ordering of the transverse
momenta in BFKL kinematics $k_{T}^{{}^{\prime}}-k_{T}\backsimeq 0$ implies the
higher order terms
$\mathcal{O}\left(k_{T}^{{}^{\prime}2}-k_{T}^{2}\right)\rightarrow 0$, thereby
ensures the higher order terms to become insignificant and can be neglected.
Thus this assumption would hold good as long as no ordering condition of the
transverse momenta in BFKL kinematics is concerned. This type of series
expansion of distribution function is well supported in the literature 12 .

Now neglecting these higher order terms we can express (16) as

$\begin{split}-x\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial
x}=&\xi(k_{T}^{2})\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial
k_{T}^{2}}+\zeta(k_{T}^{2})f\left(x,k_{T}^{2}\right)\\\
&-\Delta(k_{T}^{2})f^{2}\left(x,k_{T}^{2}\right),\end{split}$ (18)

where

$\begin{split}\xi(k_{T}^{2})=\frac{\alpha_{s}N_{c}}{\pi}k_{T}^{2}\bigg{[}&\int_{k_{T_{\min}}^{{}^{\prime}2}}^{k_{T}^{2}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{k_{T}^{{}^{\prime}2}-k_{T}^{2}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&+\int_{k_{T}^{2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{k_{T}^{{}^{\prime}2}-k_{T}^{2}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\bigg{]},\\\
\end{split}$ (19)

$\begin{split}\zeta(k_{T}^{2})&=\frac{\alpha_{s}N_{c}}{\pi}k_{T}^{2}\bigg{[}-\int_{k_{T}^{2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&+\int_{k_{T}^{2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}+\int_{k_{T_{\min}}^{{}^{\prime}2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]},\end{split}$
(20)

and $\Delta(k_{T}^{2})=\frac{18\alpha_{s}^{2}}{\pi
k_{T}^{2}R^{2}}\frac{N_{c}^{2}}{N_{c}^{2}-1}$.

Note that the integrals in (20),
$I_{1}=\int_{k_{T}^{2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}$
and
$I_{2}=\int_{k_{T}^{2}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}$
are improper since both blow-up at the lower limit $k_{T}^{2}$. To get rid of
the singularity in $I_{1}$ we have performed some angular integral
prescription. To be specific we have used the standard trigonometric integral

$\int_{0}^{2\pi}\frac{\text{d}\theta}{p+q\cos\theta}=\frac{2\pi}{\sqrt{p^{2}+q^{2}}},$
(21)

valid for $p+q>0$. Replacing the variable $k_{T}^{{}^{\prime}}\rightarrow
k_{T}^{{}^{\prime}}+k_{T}$ one can write

$\int\frac{\text{d}^{2}k_{T}^{{}^{\prime}}}{(k_{T}^{{}^{\prime}}-k_{T})^{2}}\frac{1}{(k_{T}^{{}^{\prime}}-k_{T})^{2}+k_{T}^{{}^{\prime}2}}=\int\frac{\text{d}^{2}k_{T}^{{}^{\prime}}}{k_{T}^{{}^{\prime}2}}\frac{1}{k_{T}^{{}^{\prime}2}+(k_{T}^{{}^{\prime}}+k_{T})^{2}}.$
(22)

Now using (21) in the r.h.s. of (22) we get

$\begin{split}\int\frac{\text{d}k_{T}^{{}^{\prime}}k_{T}^{{}^{\prime}}\text{d}\theta}{k_{T}^{{}^{\prime}2}(2k_{T}^{{}^{\prime}2}+k_{T}^{2}+2k_{T}^{{}^{\prime}}k_{T}\cos\theta)}&=2\pi\int\frac{\text{d}k_{T}^{{}^{\prime}}k_{T}^{{}^{\prime}}}{k_{T}^{{}^{\prime}2}\sqrt{4k_{T}^{{}^{\prime}4}+k_{T}^{4}}}\\\
&=\pi\int\frac{\text{d}k_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}\sqrt{4k_{T}^{{}^{\prime}4}+k_{T}^{4}}}.\end{split}$
(23)

Similarly using (21) in the l.h.s. of (22) we get

$\begin{split}\int\frac{\text{d}^{2}k_{T}^{{}^{\prime}}}{(k_{T}^{{}^{\prime}}-k_{T})^{2}}\frac{1}{(k_{T}^{{}^{\prime}}-k_{T})^{2}+k_{T}^{{}^{\prime}2}}&=\pi\int\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&-\pi\int\frac{\text{d}k_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}\sqrt{4k_{T}^{{}^{\prime}4}+k_{T}^{4}}}.\end{split}$
(24)

Equating (23) and (24) we obtain

$\int\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}=2\int\frac{\text{d}k_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}\sqrt{4k_{T}^{{}^{\prime}4}+k_{T}^{4}}},$
(25)

which turns out to be well-behaved for both integration limit for $I_{1}$. On
the other hand, in the limit $k_{T_{\text{min}}}^{{}^{\prime}2}\leq
k_{T}^{{}^{\prime}2}\leq k_{T}^{2}$ i.e. for not too large
$k_{T}^{{}^{\prime}2}$, the contribution from the longitudinal component to
the gluon virtuality $k^{{}^{\prime}2}$ becomes negligible which in turn
preserves the no ordering condition of transverse momentum of BFKL kinematics
i.e. $k_{T}^{{}^{\prime}2}\approx k_{T}^{2}$ very strictly. Therefore, in this
limit $k_{T_{\text{min}}}^{{}^{\prime}2}\leq k_{T}^{{}^{\prime}2}\leq
k_{T}^{2}$, it is justified to implant a factor
$\left(k_{T}^{2}/k_{T}^{{}^{\prime}2}\right)^{\lambda}$ inside the integrals
which will in fact make our calculation simpler without altering the
underlying physics. Moreover, this will fix the infrared divergence problem of
the second improper integral $I_{2}$ by allowing us to evaluate the integral
down to $k_{T_{\text{min}}}^{{}^{\prime}2}$.

Now considering these approximations in (19) and (20) and putting $I_{1}$ in
(20) we obtain

$\displaystyle\begin{split}\xi(k_{T}^{2})=&\frac{\alpha_{s}N_{c}k_{T}^{2}}{\pi}\bigg{[}\int_{k_{T_{\min}}^{{}^{\prime}2}}^{k_{T}^{2}}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{k_{T}^{{}^{\prime}2}-k_{T}^{2}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&+\int_{k_{T}^{2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{k_{T}^{{}^{\prime}2}-k_{T}^{2}}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\bigg{]}\\\
=&\frac{\alpha_{s}N_{c}k_{T}^{2}}{\pi}\frac{1}{\lambda}\left(2-k_{T}^{2\lambda}\right)\approx-\frac{\alpha_{s}N_{c}}{\pi\lambda}(k_{T}^{2})^{\lambda+1},\end{split}$
(26)

$\displaystyle\begin{split}\zeta(k_{T}^{2})=&\frac{\alpha_{s}N_{c}}{\pi}k_{T}^{2}\bigg{[}\int_{k_{T_{\text{min}}}^{{}^{\prime}2}}^{\infty}\left(\frac{k_{T}^{2}}{k_{T}^{{}^{\prime}2}}\right)^{\lambda}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\left|k_{T}^{{}^{\prime}2}-k_{T}^{2}\right|}\\\
&-\int_{k_{T_{\min}}^{{}^{\prime}2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\frac{1}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\bigg{]}\\\
=&\frac{\alpha_{s}N_{c}}{\pi}\bigg{[}\frac{k_{T}^{2\lambda}}{\lambda}-2^{2-\frac{\lambda}{2}}\lambda^{-1}k_{T}^{2\lambda}\bigg{(}1-\frac{\sqrt{k_{T}^{4}+4}}{k_{T}^{2}}\bigg{)}^{\lambda/2}\,_{2}F_{1}\\\
&-\lambda\ln\left(\frac{k_{T}^{2}}{2}+\sqrt{1+\frac{k_{T}^{4}}{4}}\right)\bigg{]}\approx\frac{\alpha_{s}N_{c}}{\pi}\left(\epsilon+\frac{(k_{T}^{2})^{\lambda}}{\lambda}\right),\end{split}$
(27)

where
$\,{}_{2}F_{1}$=$\,{}_{2}\text{F}_{1}\left(-\frac{\lambda}{2},\frac{\lambda}{2};1-\frac{\lambda}{2};\frac{k_{T}^{2}+\sqrt{k_{T}^{4}+4}}{2k_{T}^{2}}\right)$
is a standard hypergeometric function. In the above calculations we have taken
infrared cutoff $k_{T_{\text{min}}}^{{}^{\prime}2}=1\text{ GeV}^{2}$ since for
unified BFKL-DGLAP framework this provided a very consistent result towards
HERA DIS data for proton structure function $F_{2}$ 14 . In (27)
$(k_{T}^{2})^{\lambda}/\lambda$ is the only dominant term since other terms
are significantly small. The contribution from the logarithmic term in (27) is
negligible in comparison to the net contribution for all $k_{T}^{2}$. However,
phenomenological studies shows the term involving hypergeometric function in
(27) becomes irrelevant towards change in $k_{T}^{2}$ i.e. it possess a
constant value $(\approx-4.79)$ (see Fig. 3(a)). This constant contribution is
insignificant to the net contribution for $\zeta(k_{T}^{2})$ if $k_{T}^{2}$ is
enough high. But for small $k_{T}^{2}$ this contribution can not be neglected
since at this range, the $k_{T}^{2\lambda}/\lambda$ contribution itself is
small. In consequence, this constant contribution can be treated as a small
perturbation $\epsilon$ to the dominant term $k_{T}^{2\lambda}/\lambda$ . We
have performed phenomenological determination of the constant perturbation
parameter $\epsilon$ using standard non-linear regression method (see Fig.
3(b)).

Figure 3:  (a) $k_{T}^{2}$ dependence of the small perturbation
$\tilde{\epsilon}=|\epsilon|$ corresponding to the term with the
hypergeometric function in (17) (left). (b) A phenomenological calculation of
$\epsilon$ using standard non-linear regression (right). Comparison between
$\tilde{\zeta}=\frac{\pi}{\alpha_{s}N_{c}}\zeta$ (green dotted line) and
$\epsilon+2\lambda^{-1}k_{T}^{2\lambda}$ (red dotted line) is shown.

Now we are set to solve our original equation (18) 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.

Substitution of $f(x,k_{T}^{2})$ by its inverse function $\omega(x,k_{T}^{2})$
i.e. $f=\omega^{-1}$ in (18) yields

$-x\frac{\partial\omega}{\partial x}=\xi\frac{\partial\omega}{\partial
k_{T}^{2}}-\zeta\omega+\Delta,$ (28)

which is in fact a linear PDE in $\frac{\partial\omega}{\partial x}$,
$\frac{\partial\omega}{\partial k_{T}^{2}}$ and $\omega$. Now we construct a
new set of co-ordinate $\sigma\equiv\sigma(x,k_{T}^{2})$ and
$\eta\equiv\eta(x,k_{T}^{2})$ such that it transforms (28) into an ODE. To be
specific we define this transformation in such a way that it is one to one for
all $(x,k_{T}^{2})$ in some set of points D in $x$-$k_{T}^{2}$ plane. This
will allow us to solve (28) for $x$ and $k_{T}^{2}$ as functions of $\sigma$
and $\eta$. The only requirement is that we should ensure the Jacobian of the
transformation does not vanish i.e.
$\text{J}=\left|\begin{array}[]{cc}\sigma_{x}&\eta_{x}\\\
\sigma_{k_{T}^{2}}&\eta_{k_{T}^{2}}\end{array}\right|\neq 0$ in D.

Next we want to recast (28) in $(\sigma,\eta)$ plane computing the derivatives
via chain rule:

$\begin{split}&\omega_{x}=\omega_{\sigma}\sigma_{x}+\omega_{\eta}\eta_{x},\\\
&\omega_{k_{T}^{2}}=\omega_{\sigma}\sigma_{k_{T}^{2}}+\omega_{\eta}\eta_{k_{T}^{2}}.\end{split}$
(29)

Substitution of (29) into (28) yields

$-(x\sigma_{x}+\xi\sigma_{k_{T}^{2}})\omega_{\sigma}-(x\eta_{x}+\xi\eta_{k_{T}^{2}})\omega_{\eta}+\zeta\omega-\Delta=0.$
(30)

Since we want above equation to be expressed as an ODE, we require either,

$x\sigma_{x}+\xi\sigma_{k_{T}^{2}}=0\text{ or
}x\eta_{x}+\xi\eta_{k_{T}^{2}}=0.$

If we consider $x\sigma_{x}+\xi\sigma_{k_{T}^{2}}=0$, the solution $\sigma$ is
constant along the curves that satisfy

$\frac{\text{d}k_{T}^{2}}{\text{d}x}=\frac{\xi}{x}\implies\ln(xC_{1})=\int\frac{\text{d}k_{T}^{2}}{\xi}\implies
C_{1}=\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x},$ (31)

where $l=-\frac{\alpha_{s}N_{c}}{\pi\lambda}$ and $C_{1}$ is the constant of
integration. Now we can choose the new coordinates as
$\sigma=\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x}$ and $\eta=x$. Here
$\text{J}=\left|\begin{array}[]{cc}\sigma_{x}&\eta_{x}\\\
\sigma_{k_{T}^{2}}&\eta_{k_{T}^{2}}\end{array}\right|\neq 0$ as required.

We rewrite (29) as

$\displaystyle\begin{split}&\omega_{x}=\omega_{\sigma}\sigma_{x}+\omega_{\eta}\eta_{x}=-\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x^{2}}\omega_{\sigma}+\omega_{\eta}\text{,}\\\
&\omega_{k_{T}^{2}}=\omega_{\sigma}\sigma_{k_{T}^{2}}+\omega_{\eta}\eta_{k_{T}^{2}}=\frac{(k_{T}^{2})^{-1-\lambda}e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{lx}.\end{split}$
(32)

Now putting (32) in (30) we get

$\displaystyle-\eta\omega_{\eta}+\zeta\omega-\Delta=0.$ (33)

Equation (33) is an ODE and it can be solved using standard ODE solving
techniques. Now solving (33) and then transforming it to the original
coordinates $(\sigma,\eta)\rightarrow(x,k_{T}^{2})$ we get the general
solution of the KC improved MD-BFKL equation,

where
$\Gamma\left(\frac{l+n}{l\lambda}+1,-\frac{k_{T}^{-2\lambda}m}{l\lambda}\right)$
is a standard Gamma function and
$\text{G}\left(\frac{e^{-\frac{k^{-2\lambda}}{l\lambda}}}{x}\right)$ is an
arbitrary continuously differentiable function. The parameters
$m=\frac{\epsilon\alpha_{s}N_{c}}{\pi}$ and
$n=\frac{\alpha_{s}N_{c}}{\pi\lambda}$ are coming from (27). 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
(LABEL:26).

####  2.2.1 $x$ and $k_{T}^{2}$ 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_{T}^{2}$ dependence of gluon distribution setting input distribution at
some low $k_{T}^{2}$.

Let us rewrite (LABEL:26) rearranging a bit

First we will try to evaluate the functional form of the arbitrary
differentiable function G applying some initial boundary condition on it. For
$k_{T}^{2}$ evolution we set the initial distribution at $(x,k_{0}^{2})$ where
$x$ is fixed throughout the evolution. At $(x,k_{0}^{2})$ let the argument of
G be denoted by $t=\frac{e^{-\frac{k_{0}^{-2\lambda}}{l\lambda}}}{x}$.

Now writing (LABEL:27) for $(x,k_{0}^{2})$ and putting
$x=\frac{e^{-\frac{k_{0}^{-2\lambda}}{l\lambda}}}{t}$ we get

where $F_{0}$ is the unintegrated gluon distribution at $(x,k_{0}^{2})$.
Equation (LABEL:28) 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=\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x}$ in (LABEL:28) which
gives us the l.h.s. of (LABEL:27) and then we solve (LABEL:27) for
$f(x,k_{T}^{2})$. The solution for $k_{T}^{2}$ evolution turns out to be

where

$\displaystyle\begin{split}&\text{A}[x,k_{T}^{2}]=\pi
R^{2}x^{m}m^{\frac{\lambda l+l+n}{\lambda
l}}e^{\frac{m\big{(}k_{T}^{-2\lambda}+k_{0}^{-2\lambda}\big{)}}{\lambda
l}};\\\ &\text{B}[i]=18\alpha_{s}^{2}x^{m}l^{\frac{l+n}{\lambda
l}}\lambda^{\frac{l+n}{\lambda
l}}k_{0}^{2n/l}\Gamma\left.\left(\frac{l+n}{l\lambda}+1,-\frac{m(i)^{-\lambda}}{l\lambda}\right)\right|_{i=k_{T}^{2},k_{0}^{2}}.\end{split}$

Figure 4: $k_{T}^{2}$ evolution of unintegrated gluon distribution
$f(x,k_{T}^{2})$. Our results of KC improved MD-BFKL are shown for
conventional $R=5\text{ GeV}^{-1}$ (gluons are distributed throughout the
nucleus) and $R=2\text{ GeV}^{-1}$(at ”hot-spots” within proton disk).
Prediction from modified BK equation is plotted for comparison.

Figure 5: $Q^{2}$ evolution of collinear gluon distribution $xg(x,Q^{2})$. Our
results of KC improved MD-BFKL are shown for conventional $R=5\text{
GeV}^{-1}$ (gluons are distributed throughout the nucleus) and $R=2\text{
GeV}^{-1}$(at ”hot-spots” within proton disk). Theoretical prediction is
compared with global datasets NNPDF 3.1sx and CT 14.

Similarly setting the initial distribution at $(x_{0},k_{T}^{2})$, we obtain
the solution for x-evolution as follows

where, $F_{0}^{{}^{\prime}}$ is the initial gluon distribution at
$(x_{0},k_{T}^{2})$,

$\begin{split}&\Gamma_{1}=\Gamma\left(\frac{l+n}{l\lambda}+1,-\frac{k_{T}^{-2\lambda}m}{l\lambda}\right),\text{
}\\\
&\Gamma_{2}=\Gamma\left(\frac{l+n}{l\lambda}+1,m\left(-\frac{k_{T}^{-2\lambda}}{l\lambda}-\ln\frac{x}{x_{0}}\right)\right),\\\
&\text{C}[x,k_{T}^{2}]=\pi R^{2}x^{m}l^{-\frac{1}{\lambda}-\frac{n}{\lambda
l}}\lambda^{-\frac{1}{\lambda}-\frac{n}{\lambda
l}}k_{T}^{2\left(\frac{n}{l}-\lambda\left(\frac{1}{\lambda}+\frac{n}{\lambda
l}\right)\right)}\\\ &\times m^{\frac{2}{\lambda}+\frac{2n}{\lambda
l}+1}\left(-\frac{k_{T}^{-2\lambda}}{\lambda
l}-\ln\frac{x}{x_{0}}\right)^{\frac{1}{\lambda}+\frac{n}{\lambda l}},\\\
&\text{D}[x,k_{T}^{2}]=18\alpha_{s}^{2}N_{c}^{2}x_{0}^{m}l^{-\frac{n}{\lambda
l}}\lambda^{-\frac{n}{\lambda l}}e^{-\frac{mk_{T}^{-2\lambda}}{\lambda l}}\\\
&\times m^{\frac{1}{\lambda}+\frac{n}{\lambda
l}}\left(-\frac{k_{T}^{-2\lambda}}{\lambda
l}-\ln\frac{x}{x_{0}}\right)^{1/\lambda}.\end{split}$

Figure 6: $x$ evolution of unintegrated gluon distribution $f(x,k_{T}^{2})$.
Our results of KC improved MD-BFKL are shown for conventional $R=5\text{
GeV}^{-1}$ (gluons are distributed throughout the nucleus) and $R=2\text{
GeV}^{-1}$(at ”hot-spots” within proton disk). Prediction from modified BK
equation is plotted for comparison.

Figure 7: $x$ evolution of collinear gluon distribution $xg(x,k_{T}^{2})$. Our
results of KC improved MD-BFKL are shown for conventional $R=5\text{
GeV}^{-1}$ (gluons are distributed throughout the nucleus) and $R=2\text{
GeV}^{-1}$(at ”hot-spots” within proton disk). Theoretical prediction is
compared with global datasets NNPDF 3.1sx and CT 14.

We have plotted the solution for both $x$ evolution (LABEL:30) and $k_{T}^{2}$
evolution (LABEL:29) in Fig. 4 and Fig. 6. Our prediction of unintegrated
gluon distribution $f(x,k_{T}^{2})$ is contrasted with that of modified BK
equation 5 ,

$\begin{split}&-x\frac{\partial f\left(x,k_{T}^{2}\right)}{\partial x}=\\\
&\frac{\alpha_{s}N_{c}k_{T}^{2}}{\pi}\int_{k_{T_{\min}}^{{}^{\prime}2}}^{\infty}\frac{dk_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\left[\frac{f(x,k_{T}^{{}^{\prime}2})-f\left(x,k_{T}^{2}\right)}{|k_{T}^{{}^{\prime}2}-k_{T}^{2}|}+\frac{f\left(x,k_{T}^{2}\right)}{\sqrt{k_{T}^{4}+4k_{T}^{{}^{\prime}4}}}\right]\\\
&-\alpha_{s}\left(1-k_{T}^{2}\frac{\text{d}}{\text{d}k_{T}^{2}}\right)^{2}\frac{k_{T}^{2}}{R^{2}}\left[\int_{k_{T}^{2}}^{\infty}\frac{\text{d}k_{T}^{{}^{\prime}2}}{k_{T}^{{}^{\prime}2}}\ln\left(\frac{k_{T}^{{}^{\prime}2}}{k_{T}^{2}}\right)f(x,k_{T}^{2})\right]\end{split}$
(39)

which is BFKL equation supplemented by the negative nonlinear term, derived in
approximation of infinite and uniform target. In 5  the perturbative parton
saturation is studied to a vast extent, including modification of (39) in
terms of kinematic constraint, DGLAP, $\text{P}_{gg}$ splitting function and
running coupling constant, then solving the same numerically. We have also
extracted collinear gluon distribution from unintegrated gluon distribution
using the standard relation,

$xg(x,Q^{2})=\int_{0}^{Q^{2}}\frac{\text{d}k_{T}^{2}}{k_{T}^{2}}f(x,k_{T}^{2})$
(40)

sketched in Fig. 5 and Fig. 7. Our predicted collinear gluon density is
compared with that of LHAPDF global parameterization groups NNPDF 3.1sx 36
and CT 14 53 . Both of the LHAPDF datasets include HERA as well as recent LHC
data with high precision PDF sensitive measurements.

Figure 8:  (a) Comparison between linear BFKL, MD-BFKL without KC and MD-BFKL
with KC for $k_{T}^{2}$ evolution at $x=10^{-4}$ and $10^{-6}$ (b) Comparison
between linear BFKL, MD-BFKL without KC and MD-BFKL with KC for $x$ evolution
at $k_{T}^{2}=5\text{ GeV}^{2}$ and $50\text{ GeV}^{2}$. Results are shown for
$R=5\text{ GeV}^{-1}$.

Our prediction of both unintegrated and collinear gluon distribution is
obtained for two different form of shadowing: conventional $R=5\text{
GeV}^{-1}$ (order of proton radius) where gluons are spread throughout the
nucleus and extreme $R=2\text{ GeV}^{-1}$ where gluons are expected to
concentrated in hotspots within the proton disk, recalling that $\pi R^{2}$ is
the transverse area within which gluons are concentrated inside proton. From
Fig. 4-7 it is clear that shadowing correction are more prominent when gluons
are concentrated in hotspots within proton.

The $k_{T}^{2}$ (and $Q^{2}$) evolution in Fig. 4-5 is studied for the
kinematic range $1\text{ GeV}^{2}\leq k_{T}^{2}(\text{or }Q^{2})\leq
10^{3}\text{ GeV}^{2}$ corresponding to four different values of $x$ as
indicated in the figure. Our evolution for both $f(x,k_{T}^{2})$ and
$xg(x,Q^{2})$ shows a similar growth as modified BK as well as datasets
respectively for all $x$. It is also observed that the growth of
$f(x,k_{T}^{2})$ (or $xg(x,Q^{2})$) is almost linear for the entire kinematic
range of $k_{T}^{2}$ (or $Q^{2}$). This is expected since the net shadowing
term in (16) has $1/k_{T}^{2}$ dependence, which suppresses the contribution
from the shadowing term at large $k_{T}^{2}$.

The $x$ evolution of $f(x,k_{T}^{2})$ and $xg(x,Q^{2})$ is shown in Fig. 6-7
for two $k_{T}^{2}$ values viz. $5\text{ GeV}^{2}$, $50\text{ GeV}^{2}$ and
two $Q^{2}$ values viz. $35\text{ GeV}^{2}$, $100\text{ 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}\leq x\leq 10^{-2}$. We observed that the singular $x^{-\lambda}$
growth of the gluon is eventually suppressed by the net shadowing effect. KC
improved MD-BFKL seem to poses a more intense shadowing then modified BK.
However, it is hard to establish the existence of shadowing for $x\geq
10^{-3}$. The obvious distinction between the two form of shadowing $R=\text{5
GeV}^{-1}$ (conventional) and $R=2\text{ GeV}^{-1}$ (”hotspot”) is also
observed towards small-$x$ $(x\leq 10^{-3})$. Interestingly towards very
small-$x$ $(x\leq 10^{-5})$, at small $k_{T}^{2}$ (or $Q^{2}$) values (viz.
$k_{T}^{2}$=$5\text{ GeV}^{2}$, $Q^{2}$=$35\text{ 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 Fig. 8(a) we have shown a comparison between the linear BFKL equation, MD-
BFKL without kinematic constraint and MD-BFKL with kinematic constraint for
$x$ evolution. The three solutions are compared for two different values of
$k_{T}^{2}$ i.e. $5\text{ GeV}^{2}$ and $50\text{ GeV}^{2}$. Similarly in Fig.
8(b) we have shown similar comparison but for $k_{T}^{2}$ evolution for two
different values of $x$ i.e. $x=10^{-4}$ and $x=10^{-6}$. In Fig. 8 the
singular $x^{-\lambda}$ behavior of $f(x,k_{T}^{2})$ is distinct for
unshadowed linear BFKL equation. On the other hand, the deviation of 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.

####  2.2.2 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 (LABEL:26) 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_{T}^{2}$ simultaneously in $x\text{-}k_{T}^{2}$
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 9

$f(x,k_{T}^{2})=\beta\frac{x^{-\lambda}\sqrt{k_{T}^{2}}}{\sqrt{\ln\frac{1}{x}}}\exp\left(-\frac{\ln^{2}(k_{T}^{2}/k_{s}^{2})}{2\Omega\ln(1/x)}\right),$
(41)

where $\lambda=\frac{3\alpha_{s}}{\pi}28\zeta(3)$, $\zeta$ being Reimann zeta
function and $\Omega=32.1\alpha_{s}$. The nonperturbative parameter
$k_{s}^{2}=1\text{ GeV}^{2}$ and the normalization constant $\beta\sim 0.01$ 1
. Since far below saturation region both linear and nonlinear equation should
give the same solution, therefore we can take (41) 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
$\text{G}\left(\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x}\right)$
present in the general solution (LABEL:26) applying the initial distribution
(41). Let us rewrite (LABEL:27) which is the rearranged form of the general
solution (LABEL:26)

Setting initial parameters at $(x_{0},k_{T}^{2})$ we get initial distribution
(41) as

$f(x_{0},k_{T}^{2})=\beta\frac{x_{0}^{-\lambda}\sqrt{k_{T}^{2}}}{\sqrt{\ln\frac{1}{x_{0}}}}\exp\left(-\frac{\ln^{2}(k_{T}^{2}/k_{s}^{2})}{2\Omega\ln(1/x_{0})}\right).$
(43)

We denote the argument of G at $(x_{0},k_{T}^{2})$ as
$\tau=\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x_{0}}$ which implies
$k_{T}^{2}=(-l\lambda\ln(\tau x_{0}))^{-1/\lambda}$. Now writing (LABEL:33)
for $(x_{0},k_{T}^{2})$ we obtain

where, $k_{T}^{2}=(-l\lambda\ln(\tau x_{0}))^{-1/\lambda}$. Note that
(LABEL:35) is the functional form of G. We substitute
$\tau=\frac{e^{-\frac{k_{T}^{-2\lambda}}{l\lambda}}}{x}$ in (LABEL:35) which
gives us the l.h.s. of (LABEL:33) and then solve (LABEL:33) for
$f(x,k_{T}^{2})$,

Figure 9:  Three dimensional representation of (a) TMDlib HERA data fit: PB-
NLO-HERAI+II+2018 (left) (b) our solution for KC improved MD-BFKL in
$x$-$k_{T}^{2}$ phase space (right).

For simplicity let us denote r.h.s. of (LABEL:36) by $\gamma$ i.e.

$f(x,k_{T}^{2})=\gamma(x,k_{T}^{2}).$ (46)

At $x=x_{0}$ and $k_{T}^{2}=k_{0}^{2}$

$f(x_{0},k_{0}^{2})=\gamma(x_{0},k_{0}^{2}).$ (47)

Now dividing (46) by (47) we get

$f(x,k_{T}^{2})=f(x_{0},k_{0}^{2})\frac{\gamma(x,k_{T}^{2})}{\gamma(x_{0},k_{0}^{2})}.$
(48)

From (41) we have the input distribution

$f(x_{0},k_{0}^{2})=\beta\frac{x_{0}^{-\lambda}\sqrt{k_{0}^{2}}}{\sqrt{\ln\frac{1}{x_{0}}}}\exp\left(-\frac{\ln^{2}(k_{0}^{2}/k_{s}^{2})}{2\Omega\ln(1/x_{0})}\right).$
(49)

Figure 10:  Density plots showing $\lambda$ sensitivity of unintegrated gluon
distribution $f(x,k_{T}^{2})$, sketched for two canonical choice of $\lambda$
viz. $\lambda=0.4$ (left) and $\lambda=0.6$ (right).

Now substituting $f(x,k_{0}^{2})$ from (49) into (48) we obtain

$\begin{split}&f(x,k_{T}^{2})=\\\
&\frac{(\ln\frac{1}{x_{0}})^{-\frac{1}{2}}q^{n/l}k_{0}^{-2\frac{n}{l}+4\frac{l+n}{l}-2}\beta
x_{0}^{-\lambda}e^{-\frac{\ln^{2}(k_{0}^{2})}{2\Omega\ln\frac{1}{x_{0}}}}\left(k_{T}^{2}q\right)^{-\frac{l+n}{l}}\tilde{\phi}}{\left(-\frac{mq^{-\lambda}}{\lambda
l}\right)^{\frac{l+n}{\lambda
l}}\left(\tilde{\Delta}q^{n/l}\tilde{A}[k_{T}^{2}]+\chi\right)-\tilde{\Delta}k_{T}^{2n/l}\tilde{A}[q]\left(-\frac{mk_{T}^{-2\lambda}}{\lambda
l}\right)^{\frac{l+n}{\lambda l}}},\end{split}$ (50)

where

$\begin{split}\delta=&mx_{0}^{\lambda}\sqrt{\ln\frac{1}{x_{0}}}k_{0}^{2n/l-2}e^{\frac{\ln^{2}(k_{0}^{2})}{2\Omega\ln\frac{1}{x_{0}}}}\left(-\frac{mk_{0}^{-2\lambda}}{\lambda
l}\right)^{\frac{l+n}{\lambda l}},\\\
\tilde{\phi}=&\left(-\frac{mk_{0}^{-2\lambda}}{\lambda
l}\right)^{\frac{l+n}{\lambda
l}}\left(\tilde{\Delta}k_{0}^{2n/l}\tilde{A}[k_{0}^{2}]+\delta[x,k_{T}^{2}]\right)\\\
&-\tilde{\Delta}k_{0}^{2n/l}\tilde{A}[k_{0}^{2}]\left(-\frac{mk_{0}^{-2\lambda}}{\lambda
l}\right)^{\frac{l+n}{\lambda l}}.\end{split}$

Figure 11:  Density plots showing $R$ sensitivity of unintegrated gluon
distribution $f(x,k_{T}^{2})$, sketched for two form of shadowing:
conventional $R=5\text{ GeV}^{-1}$ (left) and extreme $R=2\text{ GeV}^{-1}$
(right).

Equation (50) serves as the parametric form of the solution for KC improved
MD-BFKL equation. The input distribution is inclusive in the solution (50)
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_{T}^{2}$ evolution. This parametric form of the solution actually helps us
to explore the three-dimensional insight of the gluon evolution in
$x$-$k_{T}^{2}$ phase space. However, using (50) one may also study $x$ and
$k_{T}^{2}$ evolution separately by setting $x$ fixed for $k_{T}^{2}$
evolution or $k_{T}^{2}$ fixed for $x$ evolution.

In Fig. 9 we have shown our solution of KC improved MD-BFKL equation (50) in
three dimension. The kinematic region for our study is set to be $10^{-6}\leq
x\leq 10^{-2}$ and $1\text{ GeV}^{2}\leq k_{T}^{2}\leq 10^{3}\text{ GeV}^{2}$.
In the 3D surface, the suppression in the rise of the gluon distribution
towards small $x$ due to shadowing correction is visible. However, a linear
rise of the surface in the $k_{T}^{2}$ direction can be seen, attributed to
the $1/k_{T}^{2}$ factor in the nonlinear term, which offset the effect of
shadowing at large $x$.

In Fig. 10 we have shown density plots of $f(x,k_{T}^{2})$ in
$x\text{-}k_{T}^{2}$ domain to examine the sensitivity of $f(x,k_{T}^{2})$
towards the parameter $\lambda$. Density plot allows us to visualize and
distinguish the kinematic regions with high/low $f(x,k_{T}^{2})$ in
$x\text{-}k_{T}^{2}$ plane which is more informative then any ordinary 3D
plot. In Fig. 10 the plots are sketched for two canonical choices of $\lambda$
viz. $\lambda=0.4$ and $0.6$ corresponding to two $\alpha_{s}$ values 0.15 and
0.23. Our solution seems to be very sensitive towards a small change in
$\lambda$. An apparent 50% change in $\lambda$ (0.4 to 0.6) suggests
approximately around one order of magnitude rise in gluon distribution
$f(x,k_{T}^{2})$ for an approximate limit of $x$ and $k_{T}^{2}$: $10^{-6}\leq
x\leq 10^{-5}$ and $50\text{ GeV}^{2}\leq k_{T}^{2}\leq 100\text{ GeV}^{2}$.
It is also observed that the range of high $k_{T}^{2}$ and very small-$x$ is
the high gluon distribution $f(x,k_{T}^{2})$ region where gluons are mostly
populated. In Fig. 11 we have shown $R$ sensitivity of our solution for two
choices of the shadowing parameter $R$ viz. $R=5\text{ GeV}^{-1}$ and $2\text{
GeV}^{-1}$. A satisfactory shadowing effect is observed from the comparison of
the two plots. The extreme form of shadowing $(R=2\text{ GeV}^{-1})$ is found
to suppress atleast $10\text{-}20\%$ magnitude of gluon density than the
conventional shadowing $(R=5\text{ GeV}^{-1})$ in the high $k_{T}^{2}$ and
small-$x$ region.

##  3 Equation of Critical line and prediction of saturation scale: a
differential geometric approach

Figure 12:  (a) Contour plot obtained by solving (55) (left). Solid lines show
contours of constant gradients along the curve. (b)Diagram showing critical
boundary line separating high and low gluon density region (right). Blue
dashed line is obtained from (58) while red dashed line is corresponding to
GBW model 21 .

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_{T}^{2})$
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\text{-}k_{T}^{2}$ phase space.

Although the gluon saturation can be achieved only when $Q_{s}\sim Q$ (or
$k_{T}^{2}$), 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_{T}^{2}$
and $x$ separately, the gluon distribution depends on a single dimensionless
variable $\frac{k_{T}^{2}}{Q_{s}^{2}}$ i.e.

$\phi(x,k_{T}^{2})\equiv\phi\left(\frac{k_{T}^{2}}{Q_{s}^{2}}\right).$ (51)

In recent years, this geometrical scaling property of DIS observables is
studied very extensively for various frameworks 42 ; 56 ; 57 ; 58 ; 84 . In 56
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

$-\frac{d\omega(\gamma_{c})}{d\gamma_{c}}=\frac{\omega_{s}(\gamma_{c})}{1-\gamma_{c}},$
(52)

which has been repeatedly derived in several literature 59 ; 60 ; 61 . In 42
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 5  where
a parameter $\beta$ is defined as the relative difference between the
solutions to the linear and nonlinear equation,

$\beta=\frac{|f^{\text{lin}}(x,Q_{s}^{2}(x,\beta))-f^{\text{lin}}(x,Q_{s}^{2}(x,\beta))|}{f^{\text{lin}}(x,Q_{s}^{2}(x,\beta))}.$
(53)

The crucial parameter $\beta$ actually depicts the percentage deviation that
the non linear solution shows from the linear one and it lies in the order of
$0.1\text{-}0.5$ (or $10\text{-}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 possess any
component along $f(x,k_{T}^{2})$ axis in $\mathbb{R}^{3}$. This ensures that
it does not have any direct dependence on the magnitude of gluon density
$f(x,k_{T}^{2})$, rather it depends on the rate at which $f(x,k_{T}^{2})$
changes with respect to $x$ and $k_{T}^{2}$ change. This property of gradient
actually helps us in distinguishing out the saturation region and linear
region although the distribution function $f(x,k_{T}^{2})$ has large variation
in order of magnitude for different regions in $x\text{-}k_{T}^{2}$ phase
space.

Recall that towards small-$x$, gluon evolution is suppressed due to shadowing
effect as sketch in Fig. 9 which motivates us to evaluate the gradient of
$f(x,k_{T}^{2})$ particularly along $x$ basis. For simplicity we consider an
unit vector $\vec{\nu}$ along $\tilde{\eta}\text{ }(=1/x)$ basis, we can
project the gradient $\nabla f(\tilde{\eta},k_{T}^{2})$ along $\vec{\nu}$ via
dot product $\nabla f(\tilde{\eta},k_{T}^{2})$.$\vec{\nu}$. This scalar
quantity can also be interpreted as the directional derivative along the
direction $\vec{\nu}$,

$g\equiv\nabla_{\vec{\nu}}f(\tilde{\eta},k_{T}^{2})=\nabla
f(\tilde{\eta},k_{T}^{2}).\vec{\nu}.$ (54)

Taking the Euclidean norm yields

$g=\pm|\nabla_{\vec{\nu}}f(\tilde{\eta},k_{T}^{2})|,$ (55)

where the negative (-) sign is for the descending function (or negative
slope).

We obtained a family of contours (or level curves) $\tilde{\eta}(k_{T}^{2})$
in $\tilde{\eta}\text{-}k_{T}^{2}$ plane solving (55) as shown in Fig. 12(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. 12(a) the $\tilde{\eta}\text{-}k_{T}^{2}$ 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 in 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 $\tilde{\eta}(k_{T}^{2})$ as can be seen in Fig. 12(a). In high gluon
density region, the contour curves tend to
$\tilde{\eta}(k_{T}^{2})\rightarrow\infty$ towards high $k_{T}^{2}$. 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 $\tilde{\eta}\text{-}k_{T}^{2}$ space. The level curves of the function
$g=\pm|\nabla_{\vec{\nu}}f(\tilde{\eta},k_{T}^{2})|$ are two dimensional
curves that can be obtained by setting $g=k$ where $k$ is a constant
$(k\in\mathbb{R})$. Therefore, the equations of the level curves are given by

$|\nabla_{\vec{\nu}}f(\tilde{\eta},k_{T}^{2})|=\pm k.$ (56)

Now for a known initial saturation momentum scale $Q_{s0}(1/x_{0})$ we can
predict equation of the critical line using (55) and (56). The equation of the
critical line is the equation for the level curve of the function
$g=\pm|\nabla_{\vec{\nu}}f(\tilde{\eta},k_{T}^{2})|$ that passes through the
point $(Q_{s0}(\tilde{\eta}_{0}),\tilde{\eta}_{0})$. First we find the value
of $k$ by plugging the point $(Q_{s0}(\tilde{\eta}_{0}),\tilde{\eta}_{0})$
into (56)

$k_{s0}=\pm|\nabla_{\vec{\nu}}f(Q_{s0}^{2}(\tilde{\eta}_{0}))|.$ (57)

Now the level curve passing through
$(Q_{s0}(\tilde{\eta}_{0}),\tilde{\eta}_{0})$ is obtained by setting

$|\nabla_{\vec{\nu}}f(\tilde{\eta},Q_{s}^{2})|=|\nabla_{\vec{\nu}}f(Q_{s0}^{2}(\tilde{\eta}_{0}))|,$
(58)

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 (58). In Fig. 12(b) we have
sketched a possible critical line obtained from (58) for the choice of initial
saturation scale $Q_{s0}^{2}(\eta_{0})\backsimeq 2.8\text{ GeV}^{2}$ at
$\eta_{0}=10^{6}$ (or $x_{0}=10^{-6}$) which is taken from the calculation
from the original saturation model by Golec-Biernat and Wusthoff 21 ; 5 . A
rough agreement between our prediction and that of GBW model is observed in
Fig. 12(b). However, $Q_{s}$ given by (58) is found to have weaker $x$
dependence than the one from GBW model $Q_{s}^{{}^{\prime}2}$. 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_{s}^{2}$. 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.

##  4 KC improved MD-BFKL prediction of HERA DIS data

###  4.1 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\rightarrow q\bar{q}$ process. This component from sea quark distribution
can be calculated in perturbative QCD. The relevant diagram for this QCD
process is shown in Fig. 13 and the contribution to the transverse and
longitudinal components of the structure functions can be written using the
$k_{T}$-factorization theorem 25 ; 26

$F_{\text{i}}(x,Q^{2})=\int_{x}^{1}\frac{dx^{{}^{\prime}}}{x^{{}^{\prime}}}\int\frac{dk_{T}^{2}}{k_{T}^{4}}f\left(\frac{x}{x^{{}^{\prime}}},k_{T}^{2}\right)F_{\text{i}}^{(0)}(x^{{}^{\prime}},k_{T}^{2},Q^{2}),$
(59)

where i= T, L and $\frac{x}{x^{{}^{\prime}}}$ is the fractional momentum
carried by gluon which splits into $q\bar{q}$ pair. $F_{\text{i}}^{(0)}$
includes both the quark box and cross box contribution which comes from
virtual gluon-virtual photon subprocess leading to $q\bar{q}$ production (Fig.
13 ). The gluon distribution $f(\frac{x}{x^{{}^{\prime}}},k_{T}^{2})$ in (59)
represents the sum of the gluon ladder diagrams in the lower part of the box
as shown in Fig. 13 is given by BFKL equation 2 . Here we will study the
effect of gluon shadowing using the solution $f(x,k_{T}^{2})$ of KC improved
MD-BFKL in (59).

The explicit expression for quark box contribution $F_{\text{i}}^{0}$ can be
obtained by writing four momentum in terms of the convenient light like
momenta $p$ and $q^{{}^{\prime}}\equiv q+xp$, where $x=\frac{Q^{2}}{2p.q}$ and
$Q^{2}=-q^{2}$ like as usual (see Fig. 13 ). Now we can decompose quark and
gluon momentum in terms of sudakov variables

$\begin{split}&\kappa=\alpha p-\beta q^{{}^{\prime}}+\kappa_{T},\\\
&k=ap+bq^{{}^{\prime}}+k_{T}.\end{split}$ (60)

The integration should be performed over the box diagram subject to quark
mass-shell constraints 8

$\displaystyle\begin{split}&(\alpha-x)2p.q(1-\beta)-\kappa_{T}^{2}=m_{q}^{2},\\\
&(a-\alpha)2p.q\beta-(\kappa_{T}-k_{T})^{2}=m_{q}^{2},\end{split}$ (61)

which leads to the box contribution 29

$\displaystyle\begin{split}&\tilde{F}_{T}^{(0)}(k_{T}^{2},Q^{2})=2\sum_{q}e_{q}^{2}\frac{Q^{2}}{4\pi^{2}}\int_{0}^{1}d\beta\int
d^{2}\kappa_{T}\alpha_{s}(\kappa_{T})\\\
&\times\big{\\{}[\beta^{2}+(1-\beta)^{2}]\bigg{[}\frac{\kappa_{T}^{2}}{L_{1}^{2}}-\frac{\kappa_{T}.(\kappa_{T}-k_{T})}{L_{1}L_{2}}\bigg{]}+\frac{m_{q}}{L_{1}^{2}}-\frac{m_{q}^{2}}{L_{1}L_{2}}\big{\\}},\end{split}$
(62)
$\displaystyle\begin{split}\tilde{F}_{L}^{(0)}(k_{T}^{2},Q^{2})=&2\sum_{q}e_{q}^{2}\frac{Q^{4}}{\pi^{2}}\int_{0}^{1}d\beta\int
d^{2}\kappa_{T}\alpha_{s}(\kappa_{T})\\\
&\times\beta^{2}(1-\beta)^{2}\bigg{(}\frac{1}{L_{1}^{2}}-\frac{1}{L_{1}L_{2}}\bigg{)},\end{split}$
(63)

where the denominators $L_{\text{i}}$ are

$\begin{split}&L_{1}=\kappa_{T}^{2}+\beta(1-\beta)Q^{2}+m_{q}^{2},\\\
&L_{2}=(\kappa_{T}-k_{T})^{2}+\beta(1-\beta)Q^{2}+m_{q}^{2}.\end{split}$

Figure 13:  Diagrammatic representation of the factorization formula (59)
where gluon couples to virtual photon through the (a) quark box (left) and (b)
crossed box (right) diagrams.

Note that
$\tilde{F}_{\text{i}}^{(0)}(k_{T}^{2},Q^{2})\equiv\int_{x}^{1}\frac{dx^{{}^{\prime}}}{x^{{}^{\prime}}}F_{\text{i}}^{(0)}(x^{{}^{\prime}},k_{T}^{2},Q^{2})$
i.e. the $x^{{}^{\prime}}$ integration of (59) is implicit in the
$d^{2}k_{T}^{{}^{\prime}}$ and $d\beta$ integration in (62) and (63). Now
plugging (62) and (63) in the $k_{T}$-factorization formula (59) we get

$\displaystyle\begin{split}F_{T}(x,Q^{2})=&2\sum_{q}e_{q}^{2}\frac{Q^{2}}{4\pi^{2}}\int_{k_{0}^{2}}^{\infty}\frac{dk_{T}^{2}}{k_{T}^{4}}\int_{0}^{1}d\beta\int
d^{2}\kappa_{T}\alpha_{s}(\kappa_{T})\\\
&\times\big{\\{}[\beta^{2}+(1-\beta)^{2}]\bigg{[}\frac{\kappa_{T}^{2}}{L_{1}^{2}}-\frac{\kappa_{T}.(\kappa_{T}-k_{T})}{L_{1}L_{2}}\bigg{]}\\\
&+\frac{m_{q}}{L_{1}^{2}}-\frac{m_{q}^{2}}{L_{1}L_{2}}\big{\\}}f(\frac{x}{x^{{}^{\prime}}},k_{T}^{2}),\end{split}$
(64)
$\displaystyle\begin{split}F_{L}(k_{T}^{2},Q^{2})=&2\sum_{q}e_{q}^{2}\frac{Q^{4}}{\pi^{2}}\int_{k_{0}^{2}}^{\infty}\frac{dk_{T}^{2}}{k_{T}^{4}}\int_{0}^{1}d\beta\int
d^{2}\kappa_{T}\alpha_{s}(\kappa_{T})\\\
&\times\beta^{2}(1-\beta)^{2}\bigg{(}\frac{1}{L_{1}^{2}}-\frac{1}{L_{1}L_{2}}\bigg{)}f(\frac{x}{x^{{}^{\prime}}},k_{T}^{2}).\end{split}$
(65)

Figure 14:  Prediction for proton structure function F2 obtained for two
choices of shadowing: conventional ($R=5\text{ GeV}^{-1}$) and ”hotspot”
($R=2\text{ GeV}^{-1}$). Data are taken from HERA (H1 32  and ZEUS 33 ) as
well as fixed target experiment (NMC 31  and BCDMS 34 ). Data-sets from global
parametrization groups NNPDF3.1sx 36  and PDF4LHC 15 35  is also included. The
background contribution is given by (68) with $F_{2}^{BG}(x_{0}=0.1)\approx$
0.316, 0.384, 0.391 and 0.406 corresponding to four $Q^{2}$ values viz.
$5\text{ GeV}^{2},15\text{ GeV}^{2},35\text{ GeV}^{2}$ and $45\text{
GeV}^{2}$. Separate plots of $F_{2}^{\text{npert.}}$ vs.
$F_{2}^{\text{npert.+pert.}}$ is also sketched.

Equations (64) and (65) serve as the basic tool for calculating structure
functions at small-$x$ provided that the gluon distribution $f(x,k_{T}^{2})$
is known. An analytical approach towards the calculation of $F_{2}$ and
$F_{L}$ can be found in literature 27 ; 28  for fixed coupling case using
linear BFKL equation. In the literature 9 , structure functions are calculated
taking gluon distribution $f(x,k_{T}^{2})$ from the numerical solution of
unitarized BFKL equation for running coupling consideration. It has been seen
that the analytical approach 28  considerably overestimates the actual
(numerical) solution 9 . 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 (59) $m_{q}$ denotes the quark mass and it is taken to be $m_{q}=1.28\text{
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
(62) and replacing the quark transverse momentum by
$\kappa_{T}=\kappa_{T}^{{}^{\prime}}+(1-\lambda)k_{T}$ we obtain 28

$\begin{split}&\tilde{F}_{T}^{(0)}(k_{T}^{2},Q^{2})=\\\
&\sum_{q}e_{q}^{2}\frac{\alpha_{s}}{\pi}Q^{2}k_{T}^{2}\int_{0}^{1}d\beta\int_{0}^{1}d\lambda\int_{0}^{\infty}d\kappa_{T}^{{}^{\prime}2}[\beta^{2}+(1-\beta)^{2}]\\\
&\times\lambda\frac{(2\lambda-1)\kappa_{T}^{{}^{\prime}2}+(1-\lambda)[\lambda(1-\lambda)k_{T}^{2}+\beta(1-\beta)Q^{2}]}{[\kappa_{T}^{{}^{\prime}2}+\lambda(1-\lambda)k_{T}^{2}+\beta(1-\beta)Q^{2}]^{3}}.\end{split}$
(66)

After integrating over $\kappa_{T}^{{}^{\prime}2}$ one can arrive at

$\begin{split}\tilde{F}_{T}^{(0)}(k_{T}^{2},Q^{2})=&\sum_{q}e_{q}^{2}\frac{\alpha_{s}}{\pi}Q^{2}k_{T}^{2}\int_{0}^{1}d\beta\int_{0}^{1}d\lambda\\\
&\times\frac{[\lambda^{2}+(1-\lambda)^{2}][\beta^{2}+(1-\beta)^{2}]}{\lambda(1-\lambda)k_{T}^{2}+\beta(1-\beta)Q^{2}}.\end{split}$
(67)

Equation (66) and (67) are written in terms of Feynman integral which actually
eliminates the azimuthal dependence and reduces the two fold integral
$d^{2}\kappa_{T}$ of (62) to single integral $\pi d\kappa_{T}^{{}^{\prime}2}$
. From (67) it is clear that $\tilde{F}_{T}^{(0)}(k_{T}^{2},Q^{2})$ or
$F(x^{{}^{\prime}},k_{T}^{2},Q^{2})$ possess the dimension of $k_{T}^{2}$.
Therefore, $F(x^{{}^{\prime}},$$k_{T}^{2},$ $Q^{2})$ or more conveniently
$F(x^{{}^{\prime}},k_{T}^{2},Q^{2})/k_{T}^{2}$ may be considered as the
structure function of an off mass shell gluon of approximate virtuality
$k_{T}^{2}$. In 27  differential structure functions have been studied for
fixed coupling and it is found that the ratio between longitudinal $F_{L}$ and
transverse structure function $F_{T}$ is 2:9 for fixed coupling approximation.
We have considered this ratio directly in our calculation of longitudinal
structure function. Finally, we have taken an assumption
$f(x/x^{{}^{\prime}}.k_{T}^{2})\rightarrow f(x,k_{T}^{2})$ i.e. we ignore the
$x^{{}^{\prime}}$ dependence of $f(x/x^{{}^{\prime}}.k_{T}^{2})$ which is
reasonable in LLx accuracy since

$(\ln\frac{x}{x^{{}^{\prime}}})^{n}=(\ln x)^{n}[1+O(1/\ln x)].$

The advantage of taking this assumption is that we do not have to impose the
possible constraint 8  coming from $x/x^{{}^{\prime}}<1$ on the region of
integration. In principle, the factorization formula (59) require to be run
down to $k_{T}^{2}=0$. The integral itself is infrared finite as both the
functions $f(\frac{x}{x^{{}^{\prime}}},k_{T}^{2})$ and
$F_{i}^{(0)}(x^{{}^{\prime}},k_{T}^{2},Q^{2})$ vanish at $k_{T}^{2}=0$.
However, BFKL dynamics is based on perturbative QCD which is not expected to
hold the nonperturbative small $k_{T}^{2}$ physics. On the other hand, for
small $k_{T}^{2}$ the gluon distribution vanishes linearly with the decrease
in $k_{T}^{2}$ on account of gauge invariance 9  making the contribution
small. Therefore, we have neglected this small contribution from small
$k_{T}^{2}$ region in our calculations of unintegrated gluon distribution
$f(x,k_{T}^{2})$.

Figure 15:  Prediction for proton structure function $F_{L}$ obtained for two
choices of shadowing: conventional ($R=5\text{ GeV}^{-1}$) and ”hotspot”
($R=2\text{ GeV}^{-1}$). Data are taken from HERA H1 49  and ZEUS 33 . Data-
sets from global parametrization group PDF4LHC 15 35  is also included. The
background contribution is given by (68) with $F_{L}^{BG}(x_{0}=0.1)\approx$
0.057, 0.063, 0.071 and 0.073 corresponding to four $Q^{2}$ values viz. 5
$\text{ GeV}^{2}$, 15 $\text{ GeV}^{2}$, 35 $\text{ GeV}^{2}$ and 45 $\text{
GeV}^{2}$. Separate plots of $F_{2}^{\text{npert.}}$ vs.
$F_{2}^{\text{npert.+pert.}}$ is also sketched.  Figure 16:  pQCD prediction
from KC improved MD-BFKL for the ratio $R(Q^{2})$ in the kinematic range
$2\text{ GeV}^{2}\leq Q^{2}\leq 100\text{ GeV}^{2}$. The data 49  are from H1
($\sqrt{s}=252\text{ GeV}$) and ZEUS ($\sqrt{s}=225\text{ GeV}$) experiment.
The theoretical calculation are for $\sqrt{s}=225\text{ GeV}$ analogous to
c.m.s. of ZEUS .

Before proceeding to the realistic estimation of the structure functions, we
add on the ”background” some non-BFKL contribution, $F_{\text{i}}^{\text{BG}}$
to $F_{\text{i}}$ since the above prediction of structure function is not
enough for describing DIS data 9 . This is because (64) and (65) 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_{\text{i}}^{\text{BG}}$ evolves like $x^{-0.08}$ motivated from soft
pomeron intercept $\alpha_{P}(0)=1.08$ 30 . To be precise we use

$F_{\text{i}}^{\text{BG}}(x,Q^{2})=F_{\text{i}}(x_{0},Q^{2})\left(\frac{x}{x_{0}}\right)^{-0.08}.$
(68)

An intelligent way of calculating this non-BFKL contribution is to choose
$x_{0}$ at some high $x$ and then take $F_{\text{i}}^{\text{BG}}(x_{0},Q^{2})$
from data 31  which is also listed in the figure caption. The recent HERA DIS
data taken for comparison with our result can be found in 32 ; 49  from H1
collaboration and in 33  from ZEUS collaboration. In the text 32  from H1
collaboration inclusive neutral current $e^{\pm}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\text{
GeV}^{2}$. Corresponding $F_{L}$ data of H1 is taken from 49  where
measurement are performed at c.m.s. energies $\sqrt{s}=$ 225 and 252 GeV. On
the other hand, in 33  from ZEUS collaboration reduced cross sections for $ep$
scattering for different c.m.s. energies viz. $318$, $251$ and $225\text{
GeV}$ is presented. The fixed target data from NMC 31  and BCDMS 34
collaboration which exist for $x>10^{-2}$ is also shown in the Figure 14.
Finally, we have included data from global parameterization groups viz.
NNPDF3.1sx 36  and PDF4LHC15 35  for comparison.

Figure 17:  Theoretical prediction from KC improved MD-BFKL for reduced cross
section $\sigma_{r}(x,Q^{2})$. The data are from H1 40 ; 83
($\sqrt{s}\approx$ 318, 300 and 252 $\text{GeV}^{2}$) and ZEUS 33
($\sqrt{s}\approx$ 318 $\text{GeV}^{2}$). Theoretical calculations are for
$\sqrt{s}=318\text{ GeV}^{2}$.

The $x$ dependence of the structure functions $F_{L}$ and $F_{2}$ is shown in
Fig. 14 and Fig. 15 for the four $Q^{2}$ values $5\text{ GeV}^{2}$, $15\text{
GeV}^{2}$, $35\text{ GeV}^{2}$ and $45\text{ GeV}^{2}$. The nonperturbative
contribution, $F_{2/L}^{\text{npert.}}$ (or $F_{2/L}^{\text{BG}}$) given by
(68) contrasted with total (nonpert. +pert.) contribution,
$F_{2/L}^{\text{npert.+pert.}}$. QCD prediction shows a satisfactory agreement
with DIS data for both structure functions $F_{2}$ and $F_{L}$. It is clear
that perturbative contribution is insignificant towards large $x$ ($\geq
10^{-2}$), while for small $x$ around $x\leq 10^{-3}$ this contribution
becomes visibly important. Interestingly both data and theory for $F_{2}$
structure function seem to preserve the $x^{-\lambda}$ singular behavior
rather than showing taming due to net shadowing effect in the asymptotic limit
$x\rightarrow 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\text{
GeV}^{-1}$ (gluons in hotspots) the rise of structure function becomes slower
than that of $R=5\text{ GeV}^{-1}$ ($R=R_{H}$, hadron radius) as expected. But
even the extreme form of the gluon shadowing ($R=2\text{ 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\rightarrow q\bar{q}$.
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. 15 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 $\gamma^{*}p$ cross sections of
longitudinally and transversely polarized photons $\sigma_{L}$ and
$\sigma_{T}$ as $F_{L}\propto\sigma_{L}$ and
$F_{2}\propto(\sigma_{L}+\sigma_{T})$. Since $\sigma_{L}$ and $\sigma_{T}$ are
positive, that imposes a restriction on $F_{L}$ and $F_{2}$ i.e. $0\leq
F_{L}\leq F_{2}$. To circumvent the need of a proper relationship between the
structure functions and polarized cross sections to study $\gamma^{*}p$ cross
section one can define a ratio between the structure functions or the
equivalent cross section ratio

$R=\frac{\sigma_{L}}{\sigma_{T}}=\frac{F_{L}}{F_{2}-F_{L}}.$ (69)

The advantage behind this formulation is that this is independent of any
normalization factors.

In Fig 16 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 49  where $ep$ cross section data measured for
two center of mass energies $\sqrt{s}=225$ and $252\text{ GeV}$ is enlisted.
From Fig. 16 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}\leq 5\text{
GeV})$ is also sketched in Fig. 16. 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 $\sigma_{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

$\begin{split}\frac{d^{2}\sigma^{e^{\pm}p}}{dxdQ^{2}}=\frac{2\pi\alpha^{2}}{xQ^{4}}Y_{+}\left(F_{2}-\frac{y^{2}}{Y_{+}}F_{L}\mp\frac{Y_{-}}{Y_{+}}xF_{3}\right),\end{split}$
(70)

where $Y_{\pm}=1\pm(1-y)^{2}$. The cross section for exchanged virtual boson
$(\text{Z}\text{ or }\gamma^{*})$ 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}$ $(\text{i.e. }Q^{2}\ll M_{Z}^{2}\approx 800\text{ 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 (70) can be written as

$\sigma_{r}\equiv\frac{d^{2}\sigma^{e^{\pm}p}}{dxdQ^{2}}\frac{xQ^{4}}{2\pi\alpha^{2}}\frac{1}{Y_{+}}=\left(F_{2}-\frac{y^{2}}{Y_{+}}F_{L}\right),$
(71)

where $\sigma_{r}$ is the reduced cross section. Note that (71) is symmetric
under charge exchange i.e. identical for both $e^{+}p$ and $e^{-}p$ processes
unlike (70). Additionally, it is independent of incoming electron helicity
state. Thus, at low and moderate $Q^{2}$ $(\leq m_{Z}^{2})$ 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 (69) replacing $F_{L}$ by
$F_{L}=\frac{R}{1+R}F_{2}$ which yields

$\sigma_{r}=F_{2}(x,Q^{2})\left[1-\frac{y^{2}}{Y_{+}}\frac{R}{1+R}\right].$
(72)

The $x$ dependence of $e^{\pm}p$ reduced cross section $\sigma_{r}$ calculated
for center of mass energy $\sqrt{s}=318\text{ GeV}$ is shown in Fig. 17. Our
theoretical expectation is compared with HERA H1 measurement 40 ; 83  and ZEUS
($\sqrt{s}=318\text{ GeV}$) 33 . The available H1 data for low $Q^{2}$ ($\leq
12\text{ GeV}^{2}$) measured at $\sqrt{s}=$ 318 GeV (SVX) and 300 GeV (NVX-
BST) is taken from 40 , while the same for high $Q^{2}$ ($>12\text{ GeV}^{2}$)
is taken from 83 . 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. Both theory and data agree well for our phenomenology
range $Q^{2}\leq 100\text{ GeV}^{2}$. The distinction in $\sigma_{r}$ due to
the two forms of shadowing $R=5\text{ GeV}^{-1}$ and $R=2\text{ GeV}^{-1}$ is
more prominent for smaller $Q^{2}$ values. For each $Q^{2}$, starting at some
high $x$ the reduced cross section first increases as $x\rightarrow 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\approx 0.65$ ($y=Q^{2}/sx$) and
thus characteristic turn over of cross section at $y\approx 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\rightarrow
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.

###  4.2 Virtual photon-proton cross section prediction in transition region

Traditionally 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 for exploring 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}$ 37 ; 38  in accordance with HERA data. However, for $Q^{2}\lesssim
2\text{ 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 39 . Our present study is especially focused on the
transition region $(2\text{ GeV}^{2}\leq Q^{2}\leq 10\text{ 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 40 .

Recall the neutral current $ep$ double differential cross section formula (71)
defined in the region $Q^{2}\ll M_{Z}^{2}$. 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 (71) in terms of photon-
proton reaction. In fact the structure function $F_{2}$ and $F_{L}$ in (71)
related to the longitudinally and transversely polarized photon-proton
scattering cross sections $\sigma_{L}$ and $\sigma_{T}$ by the relations

$\displaystyle F_{L}=\frac{Q^{2}}{4\pi^{2}\alpha}(1-x)\sigma_{L},$ (73)
$\displaystyle F_{2}=\frac{Q^{2}}{4\pi\alpha}(1-x)(\sigma_{L}+\sigma_{T}),$
(74)

which hold good at low $x$. Considering (73) and (74) the reduce cross section
in (71) can be written as

$\sigma_{r}=\frac{Q^{2}(1-x)}{4\pi^{2}\alpha}\sigma_{\gamma^{*}p}^{\text{eff}},$
(75)

where

$\displaystyle\sigma_{\gamma^{*}p}^{\text{eff}}=\sigma_{T}+(1-\frac{y^{2}}{Y_{+}})\sigma_{L}$
(76)

is the effective virtual photon-proton cross section. Note that the expression
for $\sigma_{\gamma^{*}p}^{\text{eff}}$ is similar to the total cross section
$\sigma_{\gamma^{*}p}^{\text{tot}}$ which is linear combination of
$\sigma_{L}$ and $\sigma_{T}$ i.e.
$\sigma_{\gamma^{*}p}^{\text{tot}}=\sigma_{L}+\sigma_{T}$. In fact the total
cross section $\sigma_{\gamma^{*}p}^{\text{tot}}$ and
$\sigma_{\gamma^{*}p}^{\text{eff}}$ can be regarded as the same quantity at
low inelasticity $y$ i.e.
$\sigma_{\gamma^{*}p}^{\text{eff}}\xrightarrow{y\rightarrow
0}\sigma_{\gamma^{*}p}^{\text{tot}}$ which differ only in the region of high
y.

Figure 18:  QCD prediction from KC improved MD-BFKL for virtual photon-proton
cross section $\sigma_{\gamma^{*}p}^{\text{eff}}$ as a function of $Q^{2}$ at
different values of $W$. The cross section for different $W$ values are
multiplied by factor multiple of 2 indicated in the figure. The included data
are from H1 40  ($\sqrt{s}=318\text{ GeV}$) and ZEUS 81  ($\sqrt{s}=300\text{
GeV}$).

Figure 18 shows the measurement of the virtual photon-proton cross section
$\sigma_{\gamma^{*}p}^{\text{eff}}$ as a function of $Q^{2}$ corresponding to
different values of $W$. The total cross section
$\sigma_{\gamma^{*}p}^{\text{tot}}$ is often expressed as a function of
$Q^{2}$ and invariant mass $W$. The standard relation between $W$, $x$ and
$Q^{2}$ is $W=\sqrt{Q^{2}(1-x)/x}$. Since $Q^{2}\approx syx$, for small-$x$ we
can have an approximate relationship between $W$ and $y$ i.e. $W^{2}\simeq sy$
which we have used in our calculations. HERA measurement for
$\sigma^{\text{eff}}_{\gamma^{*}p}$ from H1 40  ($\sqrt{s}=318\text{ GeV}$)
and ZEUS 81  ($\sqrt{s}=300\text{ GeV}$) are included for comparison with our
theoretical prediction. We have chosen $\sqrt{s}=318\text{ GeV}$ for our
calculation analogous to H1 measurement. The precision of the data is such
that their errors are hardly visible. 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. For $Q^{2}\gtrsim 3\text{
GeV}^{2}$, an excellent agreement between the theory and HERA data can be
observed for the wide range of $W$, while for $Q^{2}<3\text{ GeV}^{2}$, the
slop of the QCD prediction seems unsatisfactory. This indicates the inadequacy
of perturbative QCD at very low $Q^{2}$ ($<3\text{ GeV}^{2}$) and roughly
provides a lower bound of $Q^{2}$ to our theory.

##  5 Conclusion

In conclusion, we have presented a phenomenological study on the behavior of
unintegrated gluon distribution at small-$x$ and moderate $k_{T}^{2}$ region
particularly $10^{-6}\leq x\leq 10^{-2}$ and $2\text{ GeV}^{2}\leq
k_{T}^{2}\leq 1000\text{ 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_{T}^{2}$ dependence of unintegrated
gluon distribution function. Our prediction of gluon distribution is
contrasted with that of modified BK equation as well as global datasets NNPDF
3.1sx and CT 14. The $x$ evolution of gluon distribution shows the singular
$x^{-\lambda}$ type behavior of gluon evolution tamed by shadowing correction.
We found that for intense shadowing $(R=2\text{ GeV}^{-1})$, towards smaller
x, certainly from $x=10^{-4}$, the gluon distribution emerges from rapid BFKL
growth which indicates the dominance of gluon shadowing. While in case of
conventional shadowing $(R=5\text{ GeV}^{-1})$ an appreciable modification of
BFKL behavior can only be seen from $x=10^{-5}$. The $k_{T}^{2}$ 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_{T}^{2}$ at large values of $k_{T}^{2}$.

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 $\mathbb{R}^{3}$ into
$x\text{-}k_{T}^{2}$ phase space. We have also shown the sensitiveness of
unintegrated gluon distribution towards $R$ and $\lambda$ using density plots
in $x\text{-}k_{T}^{2}$ plane.

An important 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\otimes
F_{i}^{(0)}$. Here the unintegrated gluon distribution $f(x,k_{T}^{2})$ 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^{-\lambda}$ behavior. Even at
intense shadowing condition $(R=2\text{ GeV}^{-1})$ the $F_{2}$ structure
function is found to be suppressed only by $10\%$. In addition we have
measured $e^{\pm}p$ reduced cross section as well as equivalent $\gamma^{*}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 $\sigma_{\gamma^{*}p}^{\text{eff}}$ serves
the role of the total cross section if small inelasticity $y$ is concerned. We
have compared our predicted $\sigma_{\gamma^{*}p}^{\text{eff}}$ 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 particularly upto $Q^{2}\sim 3\text{ GeV}^{2}$ in the transition region.

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 could be a reliable framework for exploring high
energy physics over a wide range of $x$ and $k_{T}^{2}$ which is also relevant
for LHC probe and future collider phenomenology.

## Acknowledgments

We thank Prof. Johannes Bluemlein, DESY (Hamburg, Germany) and Prof. Dieter
Schildknecht, Bielefeld university (Germany) for their valuable comments. Two
of us ( P.P. and M.L.) acknowledge Department of Science and Technology (DST),
India (grant DST/INSPIRE Fellowship /2017 /IF160770) and Council of Scientific
and Industrial Research (CSIR), New Delhi (grant 09/796(0064)2016-EMR-I)
respectively for the financial assistantship.

## References

  * (1) Y.L. Dokshitzer, Zh. Eksp. Teor. Fiz 73, 1216 (1977)
  * (2) G. Altarelli, G. Parisi, Nuclear Physics B 126(2), 298 (1977)
  * (3) V. Gribov, L. Lipatov, Physics Letters B 37(1), 78 (1971)
  * (4) E.A. Kuraev, L. Lipatov, V.S. Fadin, Zhurnal Ehksperimental’noj i Teoreticheskoj Fiziki 72(2), 377 (1977)
  * (5) I. Balitsky, L. Lipatov, J. Nucl. Phys 28, 822
  * (6) M. Froissart, Physical Review 123(3), 1053 (1961)
  * (7) L.V. Gribov, E.M. Levin, M.G. Ryskin, Physics Reports 100(1-2), 1 (1983)
  * (8) A.H. Mueller, J. Qiu, Nuclear Physics B 268(2), 427 (1986)
  * (9) E.M. Levin, M.G. Ryskin, Physics Reports 189(6), 268 (1990)
  * (10) M. Devee, J. Sarma, The European Physical Journal C 74(2), 2751 (2014)
  * (11) L.V. Gribov, E.M. Levin, M.G. Ryskin, Physics Reports 100(1-2), 1 (1983)
  * (12) W. Zhu, Nuclear Physics B 551(1-2), 245 (1999)
  * (13) V. Abramovsky, Yad. Fiz. 18, 595 (1973)
  * (14) W. Zhu, K. Chai, B. He, Nuclear Physics B 427(3), 525 (1994)
  * (15) J. Ruan, Z. Shen, J. Yang, W. Zhu, Nuclear Physics B 760(1-2), 128 (2007)
  * (16) I. Balitsky, Nuclear Physics B 463(1), 99 (1996)
  * (17) Y.V. Kovchegov, Physical Review D 60(3), 034008 (1999)
  * (18) M. Kuroda, D. Schildknecht, Physical Review D 96(9), 094013 (2017)
  * (19) R. Jian-Hong, S. Zhen-Qi, Z. Wei, Chinese Physics C 32(9), 728 (2008)
  * (20) M. Ciafaloni, Nuclear Physics B 296(1), 49 (1988)
  * (21) B. Andersson, G. Gustafson, J. Samuelsson, Nuclear Physics B 467(3), 443 (1996)
  * (22) J. Kwieciński, A.D. Martin, P. Sutton, Zeitschrift für Physik C Particles and Fields 71(1), 585 (1996)
  * (23) J.R. Forshaw, P. Harriman, P. Sutton, Nuclear Physics B 416(3), 739 (1994)
  * (24) A. Askew, J. Kwiecinski, A.D. Martin, P. Sutton, Physical Review D 47(9), 3775 (1993)
  * (25) M. Kimber, J. Kwiecinski, A.D. Martin, Physics Letters B 508(1-2), 58 (2001)
  * (26) K. Kutak, A.M. Staśto, The European Physical Journal C-Particles and Fields 41(3), 343 (2005)
  * (27) V.S. Fadin, L.N. Lipatov, Physics Letters B 429(1-2), 127 (1998)
  * (28) Y.Y. Balitsky, Sov. J. Nucl. Phys. 28, 22 (1978)
  * (29) A. Askew, J. Kwiecinski, A.D. Martin, P. Sutton, Physical Review D 49(9), 4402 (1994)
  * (30) J. Kwiecinski, A.D. Martin, P. Sutton, Physical Review D 44(9), 2640 (1991)
  * (31) G. Boroun, Communications in Theoretical Physics 44(5), 891 (2005)
  * (32) A. Donnachie, P.V. Landshoff, Physics Letters B 437(3-4), 408 (1998)
  * (33) J. Soffer, O. Teryaev, Physical Review D 56(3), 1549 (1997)
  * (34) M. Devee, J. Sarma, Nuclear Physics B 885, 571 (2014)
  * (35) E. Levin, K. Tuchin, Nuclear Physics B 573(3), 833 (2000)
  * (36) E. de Oliveira, A. Martin, M. Ryskin, The European Physical Journal C 74(8), 3030 (2014)
  * (37) J. Kwiecinski, A.D. Martin, A. Stasto, Physical Review D 56(7), 3991 (1997)
  * (38) R.D. Ball, V. Bertone, S. Carrazza, L. Del Debbio, S. Forte, P. Groth-Merrild, A. Guffanti, N.P. Hartland, Z. Kassabov, J.I. Latorre, et al., The European Physical Journal C 77(10), 663 (2017)
  * (39) S. Dulat, T.J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, C.P. Yuan, Physical Review D 93(3), 033006 (2016)
  * (40) K. Golec-Biernat, M. Wüsthoff, Physical Review D 60(11), 114023 (1999)
  * (41) K. Golec-Biernat, L. Motyka, A.M. Staśto, Physical Review D 65(7), 074037 (2002)
  * (42) V. Khoze, A. Martin, M. Ryskin, W. Stirling, Physical Review D 70(7), 074013 (2004)
  * (43) A. Staśto, K. Golec-Biernat, J. Kwieciński, Physical Review Letters 86(4), 596 (2001)
  * (44) C. Marquet, L. Schoeffel, Physics Letters B 639(5), 471 (2006)
  * (45) J. Bartels, J. Blümlein, G.A. Schuler, Zeitschrift für Physik C Particles and Fields 50(1), 91 (1991)
  * (46) A.H. Mueller, D.N. Triantafyllopoulos, Nuclear Physics B 640(1-2), 331 (2002)
  * (47) S. Munier, R. Peschanski, Physical Review D 69(3), 034008 (2004)
  * (48) S. Catani, M. Ciafaloni, F. Hautmann, Nuclear Physics B 366(1), 135 (1991)
  * (49) J.C. Collins, R.K. Ellis, Nuclear Physics B 360(1), 3 (1991)
  * (50) S. Catani, M. Ciafaloni, F. Hautmann, Physics Letters B 242(1), 97 (1990)
  * (51) H1 collaboration, F. Aaron, C. Alexa, V. Andreev, S. Backovic, A. Baghdasaryan, S. Baghdasaryan, E. Barrelet, W. Bartel, O. Behrend, P. Belov, et al., The European Physical Journal C 71(3), 1579 (2011)
  * (52) ZEUS collaboration, S. Chekanov, M. Derrick, S. Magill, B. Musgrave, D. Nicholass, J. Repond, R. Yoshida, M. Mattingly, P. Antonioli, G. Bari, et al., Physics Letters B 682(1), 8 (2009)
  * (53) P. Amaudruz, M. Virchaux, R. Kaiser, M. Van der Heijden, C. Broggini, F. Sever, A. Simon, D. Nowotny, R. Schumacher, C. Mariotti, et al., Phys. Lett. B 295(CERN-PPE-92-124), 159 (1992)
  * (54) A.C. Benvenuti, D. Bollini, G. Bruni, T. Camporesi, L. Monari, F. Navarria, A. Argento, J. Cvach, W. Lohmann, L. Piemontese, et al., Physics Letters B 223(3-4), 485 (1989)
  * (55) J. Rojo, A. Accardi, R.D. Ball, A. Cooper-Sarkar, A. De Roeck, S. Farry, J. Ferrando, S. Forte, J. Gao, L. Harland-Lang, et al., Journal of Physics G: Nuclear and Particle Physics 42(10), 103103 (2015)
  * (56) J. Bartels, A. De Roeck, M. Loewe, Zeitschrift für Physik C Particles and Fields 54(4), 635 (1992)
  * (57) J. Kwiecinski, A.D. Martin, P. Sutton, Physical Review D 46(3), 921 (1992)
  * (58) H1 collaboration, V. Andreev, A. Baghdasaryan, S. Baghdasaryan, E. de Wolf, H. Jung, et al., European physical journal: C: particles and fields.-Berlin 74(4), 2814 (2014)
  * (59) A. Donnachie, P.V. Landshoff, Zeitschrift für Physik C Particles and Fields 61(1), 139 (1994)
  * (60) H1 Collaboration, F. Aaron, C. Alexa, V. Andreev, B. Antunovic, S. Aplin, A. Asmone, A. Astvatsatourov, S. Backovic, A. Baghdasaryan, E. Barrelet, et al., The European Physical Journal C 63(4), 625 (2009)
  * (61) H1 Collaboration, F.D. Aaron, C. Alexa, K. Alimujiang, V. Andreev, B. Antunovic, A. Asmone, S. Backovic, A. Baghdasaryan, E. Barrelet, W. Bartel, et al., The European Physical Journal C 64(4), 561 (2009)
  * (62) H1 Collaboration, C. Adloff, et al., The European Physical Journal C-Particles and Fields 21(1), 33 (2001)
  * (63) H1 Collaboration, C. Adloff, S. Aid, M. Anderson, V. Andreev, B. Andrieu, C. Arndt, A. Babaev, J. Bähr, J. Ban, Y. Ban, et al., Physics Letters B 393(3), 452 (1997)
  * (64) S. Donnachie, G. Dosch, P. Landshoff, O. Nachtmann, _Pomeron physics and QCD_ , vol. 19 (Cambridge University Press, 2002)
  * (65) ZEUS Collaboration, S. Chekanov et al., The European Physical Journal C - Particles and Fields 21(3), 443 (2001)