Analytical solution to DGLAP integro-differential equation in a simple toy-model with a fixed gauge coupling

We consider a simple model for QCD dynamics in which DGLAP integro-differential equation may be solved analytically. This is a gauge model which possesses dominant evolution of gauge boson (gluon) distribution and in which the gauge coupling does not run. This may be ${\cal N} =4$ supersymmetric gauge theory with softly broken supersymmetry, other finite supersymmetric gauge theory with lower level of supersymmetry, or topological Chern-Simons field theories. We maintain only one term in the splitting function of unintegrated gluon distribution and solve DGLAP analytically for this simplified splitting function. The solution is found by use of the Cauchy integral formula. The solution restricts form of the unintegrated gluon distribution as function of momentum transfer and of Bjorken $x$. Then we consider an almost realistic splitting function of unintegrated gluon distribution as an input to DGLAP equation and solve it by the same method which we have developed to solve DGLAP equation for the toy-model. We study a result obtained for the realistic gluon distribution and find a singular Bessel-like behaviour in the vicinity of the point $x=0$ and a smooth behaviour in the vicinity of the point


Introduction
DGLAP equation is a renormalization group equation (RGE) for the integrated parton distributions. It has been written initially for QED [1,2,3] in an integro-differential form. BFKL equation appears as a result of generalization of the Regge theory of scattering from quantum mechanics to QCD [4,5,6,7,8]. In Refs. [9,10] the DGLAP equation was written as a RGE for the integrated parton distributions in QCD. Dokshitzer [9] wrote this equation in an integro-differential form based on Gribov and Lipatov results in QED [1,2,3] and also Bethe-Salpeter technique used earlier in the BFKL equation was applied.
The BFKL equation is an optic theorem written down for the amplitude of scattering of two particles A(s, t) in the Regge limit. The amplitude A(s, t) may be obtained from the four-point Green function of the reggeized gluons after integrating the part of the external momenta with the impact factors. The optic theorem may be mapped to an integro-differential equation (IDE) for this four-point Green function (and for the amplitude A(s, t) in the Regge limit), in which the derivative is taken with respect to ln s. The four-point Green function depends on the variable t too. The BFKL IDE is written for unintegrated gluon distributions.
The DGLAP equation is another IDE in which the derivative is taken with respect to variable u = Q 2 /µ 2 , where Q 2 is the momentum transfer in the t-channel of the two particles in two-particles scattering process and the kernel of this IDE depends on the variables t = −Q 2 and x. The DGLAP IDE may be considered as the RGE for integrated parton distributions and is valid for large t = −Q 2 and large x in order to be in the framework of the perturbation theory. The DGLAP IDE may be written as a matrix differential equation in which the derivative is taken with respect to variable u too. This matrix differential equation is written for the Mellin moment G(N, u) of integrated gluon distribution G(x, u) and Σ(N, u) of integrated singlet distribution Σ(x, u). The procedure of the integral transformation to the Mellin moments suggests that N is a complex variable. In Refs. [11,12,13,14,15,16,17,18] an approximation of DGLAP matrix differential equation has been considered in which the Mellin moment Σ(N, u) of the integrated singlet distribution was discarded and the Mellin moment of integrated gluon distribution G(x, u) was considered only. In these articles the saddle point method was used to find an approximate solution of DGLAP equation for the Mellin moments.
The BFKL and DGLAP equations are unstable under radiative corrections in the different regimes. For example, DGLAP splitting functions are unstable at small x and BFKL kernel is unstable at large momentum transfer Q 2 . In Refs. [11,12,13,14,15,16,17,18] it has been shown that the both IDEs may be considered together on the same footing and the problem of their stability has been treated.
In Refs. [19,20] the relation between DGLAP and BFKL equations were studied from a different point of view. In N = 4 supersymmetric Yang-Mills theory due to the vanishing of β function the DGLAP splitting functions are stable even for small x and the matrix of anomalous dimensions may be obtained from the BFKL equation [20]. In N = 4 supersymmetric Yang-Mills theory the matrix of anomalous dimensions was obtained explicitly from the BFKL eigenvalues [20] without making any conjecture about the form of the splitting functions P which stand in the integral kernels in DGLAP. A possibility to obtain a matrix of anomalous dimensions from the BFKL eigenvalues in nonsupersymmetric QCD was considered in [11,12,13,14,15,16,17,18].

Integral Transforms
In this section we collect necessary formulas of various integral transforms which we will use in all the paper.

Mellin transform
We define Mellin transform as in which the arguments in the brackets on the l.h.s. stand for the transforming function f (x) and the integration variable x of this integral transformation. The inverse Mellin transformation is The position point c of the vertical line of the integration contour in the complex plane must be in the vertical strip c 1 < c < c 2 , the borders of the strip are defined by the condition that two integrals must be finite. This means that The contour in the complex plane is the vertical line with Re z = c is in the strip 0 < c < A, where A is a real and positive number, the contour must be closed to the left infinity. We may write many parameters (for example, other complex variables), − → α = (α 1 , . . . , α n ) on which the function f may depend, The integral on the r.h.s. of Eq.(1) may be seen as a sum of two integrals

Laplace transform
Representation (7) of the Mellin transformation is closely related to the Laplace transformation. We define Laplace transform of function f (x) as 2 This transformation is defined only for the functions that have restricted exponential growth a, that is f (x) < Ae ax , A is a real positive, in the right complex half-plane Re z > a. In this case the inverse transformation is where Re(z) = a + δ and δ → +0. This means that the vertical line of the integration in the complex plane passes slightly to the right of the point a. To show the compatibility explicitly, we perform subsequent transformations and obtain identity where Re z > a + δ > a. The contour is closed to the right complex infinity. We cannot close the contour to the left infinity since L[f (x), x](z) has poles in the half-plane to the left from the vertical line which crosses the real axis at the point a + δ. The inverse Laplace transformation can be checked as this is valid due to the following integral relation

Mellin moments
We define Mellin z-moment of function f (x) as z is a complex variable. To construct the inverse transformation, we need to rewrite (13) in the form of the Laplace transformation (8) and then to use (9), where we have introduced a new function F (y) ≡ f (− ln y). The Laplace transform for the function f (x) appears to be a Mellin moment for the function F (y), where a is an index of the exponential growth of the function f (x), the Mellin moment M [F (y), y](z) is defined in the same domain because the power-like restriction on its growth comes from the restrictions on f (x). In the inverse transformation the contour passes vertically in the complex plane z in the same position at Re z = a + δ as it does for the Laplace transformation (9). Under this condition the M [F (y), y](z) does not have poles in the complex half-plane to the right from this vertical line. The direct proof of the inverse transformation may be done as and the inverse transformation may be proved as where 0 < y < 1. Thus, the transformation (15) is inverse to transformation (13) under the restriction for the power-like growth (16). The Mellin moments, Laplace transform and Mellin transform posses the same equation for the inverse transformation. However, they are related by complex diffeomorphisms.

Description of Theoretical Setup
Structure functions of nucleons may be measured in deep inelastic scattering processes. They are related to integrated parton distributions which have probabilistic interpretation. There are two integro-differential equations (IDEs) for parton distributions, DGLAP equation [9,10] and BFKL equation [4,5,6,7,8], studied widely in many papers. Our paper is dedicated to a toy-model for evolution of integrated gluon distribution. We present an analytical solution to DGLAP equation in this model.

Evolution Equations
Because we need to use the evolution equations along all the paper, we briefly review the main idea of the probabilistic interpretation of them along the line of Ref. [10]. The IDE in which participates the splitting function P (x) is We calculate Mellin N -moment of both the parts of this equation and obtain the relation where we define γ(N, α(u)) as Thus, the RGE may be re-written in the form of IDE (19). A complex variable 3 N appears in the Mellin moment

Parton distributions
In the realistic QCD dynamics when there are integrated quark distributions q i (x, u) of different flavors i and there is integrated gluon distribution G(x, u), the system evolves according to IDEs given in [10], is called integrated singlet distribution and ∆ ij (x, u) are called non-singlet integrated quark distributions. Splitting functions P ab give the probability to find a parton a inside a parton b. The splitting functions may be calculated from the Lagrangian of QCD. Taking Mellin moments of both the parts of IDEs (24) and (25) we obtain matrix differential equation with the anomalous dimension matrix γ ab (N, α), where N is a complex variable which corresponds to the Mellin moments. In the present article we take into account integrated gluon distribution G(x, u) only. This approximation is known as a dominant eigenvalue of the matrix of anomalous dimensions [18,11,12,13,14,15,16,17] and may be justified in several gauge models.
The DGLAP IDE (23) may be written as a differential equation for any pair ij after taking the Mellin moments of both the parts of it in analogy with Eq. (22). As it follows from (22), the anomalous dimension γ qq (N, α) is Mellin moment of the splitting functions P qq ,

About the model and DIS processes in this model
Progress in the solution to DGLAP and BFKL equations has been achieved in N = 4 supersymmetric Yang-Mills theory [19,20]. This is due to the fact that the gauge β-function vanishes in all loops in this theory. If supersymmetry in this model is softly broken, it would not spoil the vanishing of the gauge β-function at the scale well above a gluino mass. However, the presence of a gaugino mass may make superpartners heavy while the gluons remain massless [21,22,23,24,25,26,27,28] This would mean there is no running of the coupling in the model at the scale well above the threshold of gluino mass and due to this the confinement is not possible and the existence of nuclei is doubtful in N = 4 supersymmetric Yang-Mills theory. However, the bound states of three gluinos are possible in this model, they may serve as nuclei in the analysis of DIS processes in this field theory. At the scale well below the gluino mass threshold we have a pure QCD theory without fermions and with a running gauge coupling. Massless gluons may be split into other massless gluons via the splitting functions P GG , or in a superpartners via the splitting function P Gq . The partonic model is described well by DGLAP equation. In this model there are three integrated parton distributions which are gluon distribution, gluino distribution and the corresponding scalar distribution. The corresponding solution to DGLAP IDE is a mixture of three power-like functions. However, there always is dominant contribution which has a dominant power. We treat integrated gluon distribution G(x, u) as this dominant contribution and do not take into account other two contributions from gluino distribution and scalar distribution. This is a rough approximation to DGLAP IDE of N = 4 supersymmetric Yang-Mills theory. However, it is a good model for searching analytical solution to this IDE. Such an analytical solution to DGLAP IDE is found in Section 8 and Section 9 of the present article. This approximation assumes that instead of matrix of anomalous dimensions we have only one function γ(N, α). Instead of DGLAP IDE (25) for integrated gluon distribution we consider IDE (27).
The evolution of integrated gluon distribution G(x, u) is subject to the DGLAP IDE and the evolution of the unintegrated gluon distribution is subject to the BFKL IDE. Both these IDEs must be consistent when applied to the gluon distribution which must satisfy them. The BFKL IDE is valid for each of three unintegrated distributions independently, however we consider it for unintegrated gluon distribution only because in our model we consider DGLAP IDE only for the unintegrated gluon distribution. We should consider DGLAP IDE and BFKL IDE in the kinematic region in which both equations are valid. We show in the present article that DGLAP IDE is enough to find a general form of gluon distribution in the proposed toy-model and we do not need the BFKL IDE for this purpose. Also, the model described in the previous two paragraphs possesses a property that its gauge coupling does not run. There are many gauge theories that possess such a property [21,23,25,24,29,30,34].

DGLAP equation with vanishing β-function for Integrated gluon distribution
The integrated gluon distribution G(x, u) is a dimensionless function, where u = Q 2 /µ 2 , and Q 2 is the momentum transfer and µ 2 a referential momentum transfer. It was constructed as one of the coefficient functions for the decomposition of the cross sections in DIS processes in terms of the tensor structures. In the approximation described in the previous chapters this integrated gluon distribution satisfies the DGLAP IDE, that is, We take in this section α ′ (u) = 0 and in the framework of this model we have the result for integrated gluon distribution here Q 2 = µ 2 (u = 1) is a scale which corresponds to arbitrariness in solutions to differential equations. In Refs. [31,20] it is called Q 2 0 scale. We do not write any dependence on α in the integrated gluon distributions G (N, u) , however in γ(N, α) we write it explicitly. This will be useful for further expansions in terms of α.
For the brevity, in the rest of the paper we will use the notation of Refs. [11]- [18] . From the theory of the integral transformations it follows that the small x region for the dominant PDF G(x, u) corresponds to the terms singular at the point N = 1 of the Mellin moment M [G(x, u), x](N ), see Section 9.
In this case we have a power-like dependence of PDFs on the momentum transfer [20]. Usually, there are some symmetry reasons to have the gauge coupling fixed. This happens for example in N = 4 supersymmetric Yang-Mills theory [32,33], Chern-Simons non-Abelian topological Yang-Mills theory at fixed points of the renormalization group flows [34,35], finite supersymmetric Yang-Mills theories with low level of supersymmetry [30,29], softly broken finite Yang-Mills theories [21,22,23,24,26,27,28]. In N = 4 supersymmetric Yang-Mills theory twist-two operators may be combined in representations irreducible with respect to the renormalization group with the property of multiplicative renormalization [20], and even in supersymmetric theories with the lower level of supersymmetry a dominant PDF may exist in the small x limit [36,37]. We may expect that that, if an irreducible with respect to the renormalization group multiplicatively renormalizable combination of Mellin moments of PDFs contains the moment of gluon PDF, than it is dominant in the small x limit represented by the terms singular at the point N = 1 of the complex plane of the Mellin variable. A number of these involved irreducible representations of the Mellin moments of PDFs which are dominant in the region N → 1 depends on the level of symmetry of the theory in this limit for a given theory 4 .
The solution of the DGLAP equation for the running coupling for the integrated PDF is given in Appendix A and it has been partially considered in Ref. [38] This paper is mainly dedicated to the fixed coupling so that all the comments on the case of the running coupling were put in Appendices.

DGLAP equation with vanishing β-function for Unintegrated gluon distribution
It is known that integrated gluon distribution G(x, u) is related to unintegrated gluon distribution here ϕ x, k 2 ⊥ /µ 2 is the unintegrated dominant PDF. It appears that it is always possible to construct from ϕ x, k 2 ⊥ /µ 2 a function which satisfies the same DGLAP equation (27) as well as the integrated G(x, u) dominant PDF does. In Section 5 and in Appendix B we show this statement is true in both the cases of the fixed (Section 5) and of the running gauge coupling (Appendix B). We need to consider the unintegrated PDF because the dual IDE which is called the BFKL equation is written for unintegrated PDFs [4,5,6,7,8]. We may get this dual DGLAP equation (BFKL equation) via a complex diffeomorphism from the DGLAP equation, as it has been done in Ref. [38]. This means these two IDEs, DGLAP and BFKL, should be written for the same quantities that are the unintegrated PDFs.
From Eq.(29) we conclude that their Mellin moments are related too by the same integral relation where we denoted In turn, this unintegrated gluon distribution ϕ(x, k 2 ⊥ ) solves the BFKL equation. In maximally supersymmetric Yang-Mills theory together with this function other unintegrated distributions like fermionic gluino distribution and scalar distribution exist [20]. Integrated gluon distribution is dimensionless function and unintegrated gluon distribution is dimensionful function.
From Eq.(28) we obtain This simple transformation shows that the dimensionless function Q 2 ϕ N, Q 2 /µ 2 satisfies the same DGLAP equation as Mellin moments G N, Q 2 /µ 2 of integrated gluon distribution G x, Q 2 /µ 2 do, and with the same power-like solution A new function φ N, Q 2 µ 2 may be introduced for the future use This new function φ (N, u) is Mellin N -moment of the solution to the DGLAP IDE and for this IDE the domain of u is a real nonnegative u ∈ [0, ∞[. In order to uniform notation with Appendix B dedicated to the running coupling we change the normalization of the dimensionless function φ N, Q 2 /µ 2 by a factor which is a simple constant when the coupling does not run, After this renormalization, we may show that is, the shape function φ (N, 1) of the unintegrated dominant PDF is parametrized the same way as the shape function G(N, 1) of its integrated dominant PDF is. It may be shown that a self-consistency condition should be imposed on the shape function φ 1 (N ) which may be obtained directly from the DGLAP equation in its integro-differential form. In Section 7 it is shown that such self-consistency conditions may be written for the frozen and for the running coupling. These conditions almost coincide for the cases of the running and of the fixed coupling. The self-consistency condition for the shape function in the case of the frozen coupling is applied in Sections 8 and 9. The self-consistency condition in the case of the running coupling has been obtained in Appendix C by completely the same method as it has been done in the case of the fixed coupling.

Contour of the inverse transformation from N to x
The domain of variable x of φ (x, u) should include the interval x ∈ [0, 1], otherwise the transformation to Mellin moment (37) would be impossible to define. In brief, summarizing the discussion of the previous section, if we know (37) then to recover φ (x, u) when x ∈ [0, 1] we need to make the inverse transformation (15) via Cauchy formula 5 It is supposed that Mellin moment φ(N, u) is defined in the domain Re N > a, where a is an index of the power-like growth of the function φ(x, u), In the inverse transformation (39) the contour passes vertically in the complex plane N at Re N = a + δ. Under this condition the φ(N, u) does not have poles in the complex half-plane to the right from this vertical line in the complex plane of variable N.

Method to solve the DGLAP equation analytically
In this Section we propose how DGLAP IDE may be solved without making use of the BFKL equation. This may be considered as an alternative way to the approach of Refs. [11]- [18] and to the approach of Refs. [19,20]. As we have mentioned in Introduction, the use of BFKL was a trick there to get some information about possible solution to DGLAP equation. One of the motivations for these approaches was that BFKL kernel is better known than DGLAP kernel and it was more easy to calculate the BFKL kernel than to calculate the DGLAP kernel [19,20,39] at the same loop order. DGLAP IDE (35) has a solution in the form of Eq. (34) for Mellin N -moment φ(N, u). This solution does not restrict the form of function φ 1 (N ). The reason is that when we do the integration over variable x on both sides of IDE (19), we are averaging the information about x in unintegrated gluon distribution φ(x, u). After this averaging we obtain differential equation for the Mellin moments like Eqs. (22), (28) and (36).
However, we may look at DGLAP IDE at a different angle and substitute the inverse transformation (39) in DGLAP IDE (35) for unintegrated gluon distribution φ(x, u). Such a strategy should give restrictions on function φ 1 (N ), because we use pointwise information. Indeed, by doing this we obtain The integral in the bracket may be transformed to The DGLAP IDE may be written in such a form For the future use we introduce the notation The main idea is the contour integral should be put to zero in front of each power of expansion in terms of x on the right hand side of Eq. (43) for the same contour.
The method we have proposed in this Section is based on the fact that integrals of the splitting functions in the range from 0 till x (where x is Bjorken variable) are proportional to x N where N is the complex variable of the Mellin moment φ(N, Q 2 /µ 2 ) of the unintegrated dominant PDF φ(x, Q 2 /µ 2 ). Due to cancellation of this power x N with the power x −N which stands in the inverse integral transformation, we obtain an expansion in terms of integer powers of x from which we may conclude that the coefficient in front of each integer power of x must be zero. These requirements give us a set of integrals involving the Mellin moment φ(N, Q 2 /µ 2 ) of the unintegrated dominant PDF φ(x, Q 2 /µ 2 ) which must be equal to zero simultaneously. In the next Sections 8 and 9 we have substituted the inverse Mellin moment a+i∞ a−i∞ dN x −N φ(N, u) into this DGLAP equation (35) and have obtained the equation (41) for the case of the frozen coupling, which may be treated as a self-consistency condition for the shape of the PDF. In Appendix C we simply repeat this trick for the case of the running coupling.

Solution to DGLAP equation in a simple toy-model
The IDEs of the type like Eq. (19) or in particular Eq. (35) have a probabilistic interpretation and appear in many areas of applied mathematics, mathematical biology, or stochastic processes in theoretical chemistry [40]. Some of the authors of DGLAP IDE mentioned on page 321 of textbook [41] that this equation is analogous to balance equation of various gases being in chemical equilibrium. It is not necessary that there exists a quantum field theory model for any given splitting function P (z). Quantum field theory is not the unique field of application for this IDE. The existence of a wide spectrum of applications suggests that analytical solution to such a type of IDEs should be searched. The splitting function P (z) is an input for this IDEs. In this Section we take the splitting function in the simplest form of only one term in order to show that the method we have found works for solving this IDE. Almost realistic form of the splitting function P GG (z) will be considered in the next Section.
We consider in this Section the splitting function of gluons in the form With this simple splitting function we may illustrate the main idea of the method. First, according to Eq. (44) we have We have from Eqs. (38) and (44) Thus, Eq. (43) may be rewritten in this case as from which we must conclude Eq. (49) does not restrict unintegrated gluon distribution φ(x, u) completely. Indeed, as we have explained in the previous sections, our model suggests that gluon distribution is the dominant distribution in N = 4 supersymmetric Yang-Mills theory. This is a rough approximation under which we suppose that the gauge coupling does not run and gaugino and scalar distribution are not taken into account. Coefficient β 0 is the first coefficient of the gauge β function. Since β = 0 (the coupling does not run), we have β 0 = 0. Then, Eq. (49) takes the form There are many functions satisfying this condition. For example, any expansion in powers of ln x where f k (x, u) are non-singular functions of x at x = 1, would work as gluon distribution satisfying Eq. (49). We use Eq. (49) to fix point a on the real axis in the complex plane of variable N and to find function φ 1 (N ). First, we go back to Eq. (49) and expand it in power of ln u. This expansion helps to establish the value of a, indeed, According to the theory of transformation to Mellin moment described in Section 2, all the poles should be situated to the left from the point N = a in the complex plane of variable N (⇒ −1 < a), and the contour should be closed to the negative complex infinity because x ∈ [0, 1]. There are two different possibilities to guarantee zero on the r.h.s. of Eq. (52). The first possibility is that all the poles should be of second order or higher in order to avoid contribution of residues due to Cauchy formula. This means that all the poles should be at the same point. In this Section dedicated to a simple toy-model we concentrate on this first possibility. Another possibility when residues at two different points cancel each other is considered in the next Section in which we study the solution to DGLAP by this method for almost realistic splitting function P GG (z).
We have already the pole at the point N = −1 in Eq. (49). Going along the first way described in the previous paragraph in order to solve Eq. (52) we choose that where c j are arbitrary coefficients. Another conclusion is that a is situated to the right from N = −1 on the real axis because x ∈ [0, 1] and the contour should be closed to the left, that is, −1 < a. In such a case the poles at the point N = −1 will be taken into account when we use Cauchy integral formula to calculate unintegrated gluon distribution φ(x, u). We conclude that Eq. (49) is enough to fix the contour and contains good piece of information about function φ 1 (N ). The function φ 1 (N ) could have, in fact, also terms of the form n j=1 (N − N j ) −νj with Re(N k ) ≤ −1 and at least one natural power index being ν k positive nonzero, as argued in a more general context in detail in the next Section.
As an example, we may obtain the form of unintegrated gluon distribution φ(x, u) for the simplest case when φ 1 (N ) = 1/(N + 1), that is, where I 0 is the modified Bessel function. On this side we reproduce Bessel-like behaviour obtained in Ref. [31] by summation of ladder diagrams in the pure gluonic case too. However, the Bessel-like behaviour has been obtained in Ref. [31] under some approximations for the realistic gluon splitting function P GG of Eq. (57). Our toy-model gives an exact solution for the Bessel-like behaviour with the one-term splitting function (45).
To check that the function we found possesses necessary upper bounds on its behaviour with respect to variable x, we do a simple approximation In arbitrary case we obtain This is the general solution to DGLAP IDE (35) with the splitting function (45). As we may observe, the solution is not unique. There are infinitely many constants c j which appear in this solution.
Such toy-models remain to be useful practically even nowadays because may capture in a compact expression the behaviour of a given asymptotic regime in QCD, In particular, the model (54) possesses the Bessel-like behaviour with respect to square root of the product of logarithm on the Bjorken variable and logarithm of the momentum transfer in the region of the small values of x when the main contribution comes from the gluon part of the matrix DGLAP equation. Although the computational progress of the last decades is impressive (see for example Refs. [42,43,44,45]) and the perturbative solution to the DGLAP equation is already computed up to N 2 LO for the Mellin moments of parton distribution functions with full inclusion of running coupling, the approximate solutions to the DGLAP equation corresponding to simple models still may help a lot in order to estimate physical quantities in the limits in which numerical tools and solutions show bad behaviour in the practical models like QCD.
In addition to serve as a consistency check for the numerical or analytical calculation based on a powerful software, the approximate solutions to DGLAP IDE which are presented by the models considered in this Section may be used to train neural networks [46]. Indeed, global analysis of the parton distribution functions taking into account recent data from the LHC is made by several scientific groups in the world [47,48,49,50,51]. Many PDF parameters of initial parton distribution functions may be fixed from data only because they cannot be computed from first principles. The software for the fitting of the PDF parameters and for the PDF evolution is created on the principles of neural networks [49,50] which are an efficient tool to treat a big amount of data. The forms of parton distribution functions at some scale used in such a fitting procedure tend to be some combination of Euler beta functions [52,53,54,55] which than evolve from that scale according to the DGLAP integro-differential equation.
Also, these models may be used for developing alternative analytical methods to calculate the contour integrals which appear in the inverse Mellin transformation. In Ref. [46] such contour integrals have been transformed via diffeomorphism in the complex plane of the Mellin moment variable to the contour integrals of the inverse Laplace transformation of the Jacobian of the corresponding complex map. In turn, these contour integrals of the inverse Laplace transformation may be represented in terms of the Barnes integrals by deforming the Hankel contour in the complex plane [46].

Solution to DGLAP IDE in almost realistic case
In the previous Section a toy-model has been considered. The idea was to show how the method proposed in Section 7 works. The method was aimed to solve integro-differential equations of the DGLAP type, like Eq. (19) or in particular Eq. (35). These equations have a probabilistic interpretation and due to this interpretation have many practical applications in science and technology. The toy-model was chosen to be simple, it contains one term only. However, for this toy-model we have reproduced Bessel-like behaviour of the unintegrated gluon distribution of Ref. [31] in which such a kind of behaviour has been obtained via an estimative summation of the ladder diagrams in pure gluonic QCD with the gluon splitting function P GG given in Eq. (57). This gluon splitting function P GG (x) of Eq. (57) has been calculated at the one-loop level and may be found in many textbooks.
In contrast to Ref. [31], we take the δ-function term in this splitting function equal to zero. This is because the coupling in our model does not run. This model comes from maximally supersymmetric Yang-Mills theory in which supersymmetry is softly broken. The model is described in Section 3.3. In this model we take the contribution of gluon distribution only on the r.h.s. of the DGLAP IDEs and neglect the contribution of gluino and scalar distributions. This is a rough approximation under which we suppose that the gauge coupling does not run and at the same time gaugino and scalar distribution are not taken into account. The unintegrated gluon distribution looks to be the dominant distribution in this model. This would be almost realistic model. Knowing solution in this case, we may get an impression how the gluon distribution looks in a realistic model in which all three distribution would participate.
The explicit form of the realistic gluon splitting function P GG (z) may be found in any texbook dedicated to QCD or to Quantum Field Theory in general (for example in Ref. [56], page 236, Eq. (8.5.42)), or in the original paper [10], and it takes the form in which β 0 is the one-loop coefficient of the gauge β-function.
We have to put β 0 = 0 because the coupling does not run in the case that we consider in this paper. This point requires a special comment. In N = 4 supersymmetric Yang-Mills theory the coupling does not run to all the loops. However, Eq. (57) is just a leading-order contribution to the splitting function P GG (z). In the original papers of [10,9] the splitting function P GG (z) corresponds to the kernel of Bethe-Salpeter equation [9]. We do not consider higher-order corrections to the splitting function P GG (z) in the present paper. Thus, the solution to the DGLAP equation with the splitting function (57) is the solution but only at the leading order. Its order is determined by the order of the splitting function. We do not consider other splitting functions due to the reasons that we have explained in the previous Sections. The gluon distribution dominates in the small x limit in QCD and in the conformal gauge theory like N = 4 supersymmetric Yang-Mills theory.
Altarelli and Parisi in Ref. [10] have shown that the approach based on the operator product expansion used in the Nobel prize paper [57] admits a probabilistic interpretation in terms of the splitting functions (57). It was found in Ref. [10] that these splitting functions are consistent with the anomalous dimensions of the twist two operators calculated in [57]. Similar splitting functions appeared in the approach of Refs. [1,2,9] based on the Bethe-Salpeter equation imposed on the contributing family of Feynman diagrams.
The coefficient 2C 2 (G) in the expression for the splitting function (57) is actually 2N for the gauge group SU (N ) [10]. For QCD, for example, we consider the group SU (3). Thus it is a universal coefficient based on the gauge group contribution, it does not depend on the representation of the quark fields. However, the coefficient β 0 is very sensitive to the representation of the matter fields. In QCD this coefficient is responsible for the phenomenon of the asymptotic freedom [57].
The solution to the DGLAP IDE for the Mellin moment of the dominant parton distribution is given in Appendices A and B. At the leading order of the perturbation theory for the case of the running coupling the solution to the DGLAP IDE may be represented in the same form of the contour integral (54) which we obtained for the case of the fixed coupling. The only difference with the fixed coupling case is that instead of the power function of u in the integrand of (54) another dependence on the momentum transfer u will stand. At higher orders of the perturbation theory dependence of the integrand on the momentum transfer u may be more complicate.
To calculate T (N, x, α) of Eq. (44) for this model, we need to take into account that where C is Euler-Mascheroni constant. Integral (60) comes from the first term in the gluon splitting function (57) which is defined as This means that This integral generates harmonic numbers and generalizes them to the complex argument z, Here we use the well-known binomial expansion for an arbitrary complex power z and x ∈ [0, 1].
in which (a) k = Γ(a + k)/Γ(a) stands for Pochhammer symbol. This formula may be derived by using Mellin-Barnes transformation [58]. In particular case, when x = 0 we obtain for integral (63) Also, another representation of Euler digamma function necessary for future use is Taking into account that ψ(n + 1) = H n − C, where n is a natural number, integral (65) may be considered as an analytic continuation of harmonic numbers to the complex plane z. In such a case integral (65) is an analytic continuation of Euler integral According to Eq.(44) and Eqs. (58), (59) and (60) we have for T (N, x, α) with P GG (57) As we have mentioned, this Section is based on generalization of the solution for the toy-model considered in the previous Section. Thus, we should write for the anomalous dimension We note, that in this model the normalization condition cannot be maintained due to the pole in the complex plane at the point N = 1. In view of Eq.(66) we may re-write Eq.(70) Thus, in analogy to the toy-model of the previous Section Eq. (43) may be rewritten in this case as We have obtained that some infinite series of the integer powers of x must be zero. This means that the coefficient in front of each power is zero, that is, the following identity must be fulfilled In analogy to Eq. (52) of the toy-model we obtain in which ν j ∈ N ∪ {0} are arbitrary natural numbers or zero, N j belong to a set of arbitrary complex numbers such that Re N j < a, and at least one of the numbers ν j should be nonzero. The requirement Re N j < a guarantees that all the poles of function φ 1 (N ) appear to the left from the vertical line of the contour in the complex plane N. To fulfill Eq. (74) by the terms of Eq. (76) we need to require that a > 1.
To prove that a term like (76) gives a solution to Eq. (75), we consider a simplified form of φ 1 (N ) where λ ∈ C is an arbitrary complex number such that Re λ < a. We may consider a term According to the theory of transformation to Mellin moment described in Section 2, all the poles should be situated to the left from the point N = a in the complex plane of variable N and the contour should be closed to the negative complex infinity because x ∈ [0, 1]. If a > 1 than Let us consider another combination, and we obtain again The terms of second degree or higher do not contribute into residue calculus due to Cauchy formula and the first term does not contribute due to Eq. (79). The third type of terms, which we consider in this proof, is Such a representation means that At the end of this proof, we observe that any term of type like in Eq (76) may be decomposed in a finite sum of terms (77) or their natural powers. Formulas (79), (81) and (83) show that the term (77) is a solution to Eq. (75) if a > 1. Thus, any linear combination of the terms like (76) can be used for function φ 1 (N ). To show how the residue calculus works for this solution, we take again the simplest case which has been used in the previous Section for the toy-model and has appeared to be successful in reproducing the Bessel-like behaviour of unintegrated gluon distribution φ(x, u) reviewed in Re. [31].
The result of calculation for the first two orders of expansion in terms of powers The integrals may be taken by Cauchy formula in each power of ln u. The integral in front of the first power of α 2π ln u is Here we take into account the integrals and the integral a+i∞ a−i∞ As we may see in Eq.(85) there are singularities at the points x = 0 and x = 1 at the first order of the expansion in terms of (α/2π) ln u. We sum the leading terms of these singularities and show that the singularity at x = 0 survives while the singularity at x = 1 disappears. First, we treat the singularity at the point x = 0. It is produced by the residue at N = 1. The most singular contribution is produced by the natural powers of 1/(N − 1) in each term of the expansion in Eq.(85) because in addition to 1/x we will obtain factor ln (1/x) in the maximal power. Thus, in the vicinity of the point x = 0 we may write This equation gives by itself an upper bound on unintegrated gluon distribution φ(x, u) in the vicinity of the point x = 0. The upper bound is a singular function at the limit x → 0. To be sure that φ(x, u) is a singular function we need to consider a lower bound for it in the vicinity of the point x = 0, where, when using the asymptotic behavior of the modified Bessel function I 0 (z), we have and z ≡ ln(1/x) and K = (4C 2 (G)α/π) ln u.
We conclude that a lower bound for the unintegrated gluon distribution in the vicinity of the point x = 0 is determined by the modified Bessel function and it is singular in the small x region. On the contrary, the singularity at the point x = 1 disappears. We may conclude for the considerations presented in the previous paragraphs of this Section that the most singular contribution is the biggest power of ln (1 − x). This may come only from powers of H(N ) = ψ(N +1)+C Harmonic number function in Eq. (85). If we consider integral we conclude by considering carefully the singularity structure of the Harmonic number function H(N ) = ψ(N + 1) + C in the complex plane, applying repeatedly the Cauchy theorem, and then using the identities that in the vicinity of the point x = 1 the result for integral (94) has the following asymptotic behaviour The same is true for the higher power of the ψ(N + 1) function in the integrand of Eq. (85). Thus, at the vicinity of the point x = 1 we may write We observe that the highest singularities at the point x = 1 disappear after summing the leading singularities up. This is in agreement with Eq.(74). Indeed, φ(1, u) looks like Eq.(74) without the denominator in the integrand.

Conclusion
In the present article we have found a way to solve DGLAP integro-differential equation analytically. The method we propose is simple and is based on the fact that integrals of the splitting functions in the range from 0 till x (where x is Bjorken variable) are proportional to x N where N is the complex variable of the Mellin moment φ(N, Q 2 /µ 2 ) of the unintegrated gluon distribution φ(x, Q 2 /µ 2 ), cf. Eq. (41). Due to cancellation of this power x N with the power x −N which stands in the inverse integral transformation, cf. Eqs. (42)- (43), we obtain an expansion in terms of integer powers of x from which we may conclude that the coefficient in front of each integer power of x must be zero. These requirements give us a set of integrals involving Mellin moment φ(N, Q 2 /µ 2 ) of unintegrated gluon distribution φ(x, Q 2 /µ 2 ) which must be equal to zero simultaneously, cf. Eqs. (49) and (74). We have found a way to solve these integral restrictions analytically by making use of Cauchy formula. The method we have found may have a wide spectrum of applications in science and technology.
We have considered a simple toy-model of DIS processes and found an analytical solution for the DGLAP equation in this toy-model. A simplified splitting function (45) was used as an input. The Mellin moment φ(N, Q 2 /µ 2 ) of unintegrated gluon distribution φ(x, Q 2 /µ 2 ) appears to be a linear combination of the chosen terms. The infinite set of constants c j which are coefficients in front of these chosen terms remains unfixed in this toy-model. The solution is parametrized by them. It could be that they are fixed if we consider DGLAP IDE together with BFKL IDE. However, we have shown in this article that the corresponding DGLAP IDE by itself contains enough information to represent the chosen unintegrated gluon distribution φ(x, Q 2 /µ 2 ) in this toy-model in the form of expansion in terms of ln Q 2 /µ 2 and ln 1/x shown in Eq. (56).
When we choose only one simplest term from all the possible terms, we obtain a Bessel-like behaviour for unintegrated gluon distribution φ(x, Q 2 /µ 2 ). Such a behaviour of φ(x, Q 2 /µ 2 ) has been obtained in Ref. [31] by summing ladder diagrams in an estimative way for the realistic splitting function P GG (z) in a pure gluonic Chromodynamics. We have shown that such a behaviour corresponds to the selection of this simplest term from all the possible terms for the Mellin moment φ(N, Q 2 /µ 2 ) of φ(x, Q 2 /µ 2 ) in our toy-model with the simplified splitting function.
Situation becomes more complicated for the realistic one-loop splitting function P GG (z). The number of the possible terms for the Mellin moment φ(N, Q 2 /µ 2 ) is infinite too, however more rich structure of the splitting function produces more complicate anomalous dimension for unintegrated gluon distribution. As the result, the distribution φ(x, Q 2 /µ 2 ) looks more complicated than for the toy-model. This happens even in the case when the same simplest term like in the toy model is selected of all the possible terms for Mellin moment φ(N, Q 2 /µ 2 ). Making complex integrals by use of Cauchy formula for the selected simple term of the Mellin moment φ(N, Q 2 /µ 2 ), we obtain distribution φ(x, Q 2 /µ 2 ) as an expansion in powers of α ln Q 2 /µ 2 .
The summation of this expansion in powers of α ln Q 2 /µ 2 looks difficult in this realistic case. However, the second term of the expansion shows singularities at the points x = 0 and x = 1 whose origin in the complex plane of variable N may be detected and the corresponding terms responsible for these singularities may be analysed. These singular terms at the points x = 0 and x = 1 may be summed up in all the orders of the expansion in powers of α ln Q 2 /µ 2 . After summing up these singularities at the point x = 1, they disappear and the behaviour of unintegrated gluon distribution φ(x, Q 2 /µ 2 ) becomes smooth with respect to variable x in the vicinity of the point x = 1. However, the sum of the singular terms at the point x = 0 taken to all orders of α ln Q 2 /µ 2 remains singular with respect to x at the point x = 0. The result of summation shows the Bessel-like behaviour in the vicinity of x = 0 which is similar to the behavior of unintegrated gluon distribution obtained in Ref. [31] by summing ladder diagrams or by calculating integrals via saddle-point method.
We found in this paper a large set of solutions to the DGLAP equation without using any other information from any additional equation. In particular, we did not use any information from BFKL equation. This may be considered as an alternative way to the approach of Refs. [11]- [18] where BFKL IDE has been widely used. We have shown that this integro-differential equation has infinitely many solutions for any given kernel P (z) by itself if we do not provide any boundary condition for unknown parton distributions.

A Running coupling case
In the case when the gauge coupling runs, that is the case of QCD, the first order differential DGLAP equation (27) in the small x limit for the Mellin moment G (N, u) of the dominant PDF where the coupling α has been chosen as a variable in a usual way instead of the scale u for this equation, β(α) = u dα(u)/du. Then, Such a change of variable requires that α(u) is a monotonic function. This is true in the perturbation high energy QCD [59,60]. Here G (N, u(α 0 )) = G (N, 1) is a Mellin moment of the shape function at the scale Q 2 = µ 2 , that is, at u = 1, α 0 = α(1). According to our notation, G (N, 1) appears to be the Mellin moment of G(x, 1). This function G(x, 1) should be parametrized. We mentioned in the Introduction of Ref. [38] various known parametrizations of the PDF shapes at some fixed momentum transfers for the case of QCD.
We may do the same comments, that we have done in Appendix A for the solution to the differential equation for the Mellin moment G(N, u) of the integrated dominant PDF, in the QCD case, that is, in the case when the coupling runs. Namely, the change of variables from the momentum transfer u to the coupling α supposes one-to-one correspondence between u and α. This happens at least in the penetrative high energy QCD [59,60]. Here φ (N, u(α 0 )) = φ (N, 1) ≡ φ 1 (N ) is the Mellin moment of the shape function at the scale Q 2 = µ 2 , that is, at u = 1, α 0 = α(1). The moment φ(N, 1) of the shape function φ(x, 1) for the unintegrated dominant PDF may be obtained from the solution (101) to the integrated dominant PDF G(N, u) and the parametrization of the shape function G(x, 1). We have written in Ref. [38] about the parameterizations which are frequently used for the shape function. It may be proven from Eqs. (100) and (29)  This means, the shape function φ (N, 1) of the unintegrated dominant PDF coincides with the shape function G(N, 1) of the integrated dominant PDF, and they are parametrized identically.
C Self-consistent shape function for the running coupling The DGLAP IDE (103) has a solution in the form of Eq. (108) for the Mellin N -moment of the unintegrated dominant PDF φ(N, u). This solution does not restrict the form of the function φ 1 (N ). The reason is that when we do the integration over variable x on both sides of IDE (103), we are averaging the information about x in the unintegrated dominant PDF φ(x, u). After this averaging we obtain a differential equation for the Mellin moments like Eqs. (107) and (104). However, as in the case of fixed α we may look at DGLAP IDE at a different angle and substitute the inverse transformation (15) in DGLAP IDE (103) for the unintegrated dominant PDF φ(x, u). Such a strategy should give restrictions on the function φ 1 (N ), because we use pointwise information. Indeed, by doing this we obtain The main idea to get self-consistency condition is the contour integral should be put to zero in front of each power of expansion in terms of x on the right hand side of Eq. (111) for the same contour.