Numerical solution of DGLAP equations using Laguerre polynomials expansion and Monte Carlo method

We investigate the numerical solutions of the DGLAP evolution equations at the LO and NLO approximations, using the Laguerre polynomials expansion. The theoretical framework is based on Furmanski et al.’s articles. What makes the content of this paper different from other works, is that all calculations in the whole stages to extract the evolved parton distributions, are done numerically. The employed techniques to do the numerical solutions, based on Monte Carlo method, has this feature that all the results are obtained in a proper wall clock time by computer. The algorithms are implemented in FORTRAN and the employed coding ideas can be used in other numerical computations as well. Our results for the evolved parton densities are in good agreement with some phenomenological models. They also indicate better behavior with respect to the results of similar numerical calculations.

There exist several analytical and numerical methods to solve DGLAP evolution equations. What we present in this article, is a solution which is based entirely on numerical analyses of these equations, forming a series of Laguerre polynomials with respect to the y = ln(1/x) variable, where x is the fraction of proton momentum carried by the parton, The Laguerre series converges very quickly and it can easily be truncated with a reasonable precision.
In this article, we assume that reader is familiar with the required relations and the theoretical frameworks based on which the DGLAP evolution equations are working. So, at different sections of this article, we are mostly focused on presenting numerical investigations which would finally yield the evolved parton densities at energy scale Q 2 . In each section of the paper, we not only introduce the required theoretical expressions but also explain how to use them in practice to do our numerical calculations. The Monte Carlo algorithms which we construct are such as to make the wall clock time by computer in a proper time. The numerical patterns which we develop in this paper can also be used for other numerical investigations.
The organization of this paper is as following. In section "A short overview of the theoretical framework" we give a short overview on the evolution of parton densities, using DGLAP equations. The theoretical framework is based on using Laguerre polynomial expansions. In section "Basic tools, Monte Carlo solutions" the general structure of Monte Carlo algorithm is introduced. The required functions and subroutines are also introduced there. They can be requested via E-mail, a.mirjalili@yazd.ac.ir, from the authors. Then we use them to build the related Monte Carlo algorithm to get the numerical solutions for DGLAP evolution equations. The results are presented for the evolved parton densities at the end of section "Basic tools, Monte Carlo solutions" and also in sections "The programs and the results" and "Conclusion", based on the numerically Monte Carlo algorithm. The results are in good agreement with the results of CTEQ and GRV parameterization groups. Finally, we give our conclusion in section 6.

A short overview of the theoretical framework
In high energy physics, the parton densities at the Q 2 scale can be obtained, employing the DGLAP evolution equations. These equations can be used to describe Bjorken violation in deep inelastic scattering (DIS). There are many different ways to solve DGLAP equations numerically. One of them is to use the Laguerre polynomial expansion which we employ to get the solution of non-singlet and singlet sectors of parton densities. To numerically achieve the evolved parton densities, we initially need to reach high levels of precision. The Laguerre polynomial expansion are rapidly convergent for medium values of x and at all energy scales. At very lowx, say x < 0.001, these polynomials are numerically instable due to rapid rising of the splitting moments.
Following, we have provided the required definitions and conventions for the above mentioned numerical calculations:

Running coupling constant
The Q 2 dependence of the strong coupling constant α(Q 2 ), considering the renormalization group equation, is given by: where in which β 0 and β 1 are universal scheme independent coefficients and are given by: where Nand n f refers respectively to the number of quark colors and flavors. The solution of β-QCD evolution equation, Eq. (1), at the nextto-leading order (NLO) approximation is given by: The cutoff parameter Λ, is determined by fitting the experimental data which at the NLO approximation is lower than 250 MeV.

DGLAP evolution equations
Considering the contribution of quark-antiquark pair in evolution of quark densities, one would find out that for each quark flavor i with i = 1…2n f summing over all quark and antiquark flavors, using the notations of Furmanski and Petroznio (1982a, b) we would have Following that by defining and the new combination of splitting functions by Ellis et al. (1996): the evolved DGLAP equations will be appeared in the following forms: (1) T R n f , (1) where ⊗ symbol is indicating the following convolution integral: Similar expansion like Eq. (7) exists for the different elements of the related matrix of splitting function. The advantage of using Eqs. (8-10) are that we are able to extract the sea quark densities at any energy scale Q 2 separately for each quark flavor rather than to get an average quantity for sea quark densities.
Equations (8-10) can be written in terms of the new variable t = − 2 β 0 ln α(Q 2 ) α(Q 2 0 ) so as: where Solutions of Eqs. (10-12) will lead us to the evolved valence, sea and gluon densities at different energy scales.

Evolution operators
Defining q(x) ≡ q(t = 0, x) as parton density at initial energy scale Q 0 , the evolved valence density will be obtained by For the χ i function we will have and for the gluon and singlet distributions, we will have: In Eq. (19), the first term on the right hand side is the evolution operator which for the singlet and gluon densities has the following matrix form: Substituting Eqs. (17)(18)(19) in Eqs. (12-14) will lead us to: We should note that evolution operators satisfy the following initial conditions where I In Eq. (23) refers to the unit matrix with dimension 2.

Laguerre expansion
The Laguerre polynomials can be represented by an alternative form as Arfken and Weber (2005): These polynomials have the following properties: 1. The generating function of these polynomials is indicated by: 2. They satisfy the following recursive relation: 3. These polynomials possess a closure property under the convolution integral: 4. They also satisfy the following orthonormal condition, using the weight function e −y : ∞ 0 dye −y L m (y)L n (y) = δ m,n .
Since the polynomials form a complete set, any function can be expanded in terms of them: Following that we will have: Based on the completeness property, any two arbitrary functions A(y) and B(y) can be expanded in terms of Laguerre polynomials: Assuming C(y) = A(y) ⊗ B(y) = ∞ n=0 C n L n (y) and using the closure property, given by Eq. (27), the following relations can be obtained between the expansion coefficients where If we wish to use the Laguerre polynomials to get the evolved parton densities, we should change the interval of the y variable from (0, ∞) to (0, 1). Therefore, we need to change the variables as in the following Following that the orthonormal condition, given by Eq. (30), can be written as Now the expansion of an arbitrary function F(x) would take the form so as for expansion coefficient, F n , we can write Now we are equipped with the required relations to extract the evolution operators for parton densities in terms of the Laguerre polynomials which will be done in next subsection.

Evolution operator for the non-singlet density: LO approximation
At the leading order (LO) approximation, Eq. (21) for the non-singlet evolution operator can be written as Substituting Eqs. (32,33) in Eq. (38), we arrive at By using the initial condition, given by Eq. (23), for non-singlet sector, the general solution is where By substituting Eq. (42) in Eq. (41), the evolution operator for the non-singlet sector at the LO approximation is determined and we can obtain parton densities at energy scale Q 2 based on Eq. (17).

Evolution operator for the non-singlet density: NLO approximation
We intend now to obtain the solution of the following differential equation where Ris determined by Eq. (15). We can write the following Laguerre expansion for R − Similar Laguerre expansions exist for E − (t, x) which by substituting in Eq. (43), using Eqs. (32, 33) we will arrive at x . Ghasempour Nesheli et al. SpringerPlus (2016) 5:1672 where Now at the NLO approximation for E n (t), we can write To simplify the calculation we present AE (1) n (t) by S n (t). Therefore Eq. (45) can be written as Using Eq. (40) and the initial condition, given by Eq. (23), we will get In summary, the E n (t) term up to NLO approximation would have the following form where As in the LO approximation, by substituting Eq. (50) in the related expansion for E − (t, x) , the evolution operator for the non-singlet sector at the NLO approximation is determined and finally the parton densities at any energy scale can be obtained by using Eq. (17).

Evolution operator for the singlet and gluon densities
This subsection contains two parts. At first, the solutions for the singlet sector, q (+) , and gluon densities are introduced. At the next step, using the solution for singlet sector and the χ i distribution, it is possible to get the solution for q (+) i which is defined by Eq. (5). Then, by accessing the valence distribution from the previous subsection and the q (+) i distribution, sea quark distributions for individual flavors will be obtained. More details of these calculations are as followed.
In order to be able to extract sea quark densities, at first Eq. (9) should be solved and in terms of evolution operator, we will have [see Eq. (21)]: The solution of this equation at the LO approximation is like the non-singlet case. At the NLO approximation, we should just do the following replacements with respect to the non-singlet case: At the next step, Eq. (14) for singlet and gluon densities should be solved which at LO and NLO approximations could be done as follows:

LO approximation
At LO approximation, for related evolution operator, we have [see Eq. (22)] The matrix evolution operator in Eq. (54) can be expanded in terms of Laguerre polynomials, so as The general solution for elements of the matrix evolution, by analytical consideration of moments, is as follows: where e 1 , e 2 are projection matrix operators with the following properties λ 1 (s) and λ 2 (s) are eigenvalues of the P (0) (s) matrix which are given by and due to momentum conservation at each vertex of parton splitting, we will have λ 2 (1) = 0 (Furmanski and Petroznio 1982a, b).
Considering the Laguerre expansion of evolution operator and also the splitting function in the Matrix form and Eq. (22) which connect these two expansions to each other, we will arrive at: In Eq. (59), the A (k) n and B (k) n coefficients are two dimensional matrix which are obtained from the following recurrence relations (Furmanski and Petroznio 1982a, b) (53) n in Eq. (59) and the result of the related Laguerre expansion for the matrix evolution operator, the operator is obtained and we are able to use Eq. (19) to evolve singlet and gluon densities at higher energy scales. Equations (19,20) at the LO approximation can be represented by: Using Eq. (65), we can obtain gluon and singlet distribution, q (+) , and then using the evolved valence quark (non-singlet) and χ i distributions, Eqs. (17, 18), the sea distributions at the LO approximation will be obtained [see Eqs. (5, 6)]

NLO approximation
The evolution operator in Eq. (22) at the NLO approximation is written as: The LO contribution, E (0) (t, x), is given by Eq. (59) and for the NLO contribution we should use the following relations (Furmanski and Petroznio 1982a, b) e 1 ≡ e 1 (1), e 2 ≡ e 2 (1), where Ẽ (1) n can be obtained from E (0) n , using the following integral: where R j is the expansion coefficients of R(x) in terms of Laguerre polynomials so as The sum in Eq. (69) should be used, considering the condition which is given by Dirac Delta function, that is n = i + j + k.
As before by accessing as well to the non-singlet and χ i densities, using Eqs. (19, 50, 59, 68) the sea quark densities at the NLO approximation can be extracted Now equipped by the required theoretical framework, based on Laguerre polynomial expansion, we will be able to get numerical results for the parton densities at any energy scale.

Basic tools, Monte Carlo solutions
Here, we fully describe the numerical solutions of DGLAP evolution equations, using the FORTRAN programing language. At first, we introduce the general representations of programs, functions and subroutines which we are using in all our FORTRAN codes. The programs are divided into two parts, including the LO and NLO approximations; and each part contains a non-singlet and a singlet section. We then present the obtained numerical results by depicting the related parton densities at Q 2 = 4, 50 and 200 GeV 2 . We also compare our results with the related results from CETQ and GRSV phenomenological groups. Comparisons with other numerical results have also been carried out.

Functions, subroutines and main programs
We write the required codes in FORTRAN 90 language to solve numerically the DGLAP evolution equations. The basic method to solve the integrals is to use Monte Carlo simulations. The only generic subroutine used is Ran3 which generates random numbers. All the other subroutines and function are written by us. We first introduce the functions and subroutines written by us and we will then illustrate the compiled programs in different sections.

Run3(idum) function
This function generates random numbers with uniform distribution between 0 and 1 based on Park-Miller method and Knuth suggested corrections. Each negative integer number can be considered as the input idum. The generation of this input should not be changed when we call them subsequently. The order of the period interval of this generation is 10 8 (Press et al. 1996). The random numbers are generated in the interval [a, b], using Ran3 function base on the following formula Using Eq. (70) repeatedly, we can see that the generated numbers have uniform distribution and therefore we can use them to perform Monte Carlo integrations.
Monte Carlo integration This method of integration is used to obtain the definite integrals which are based on generating random numbers. There are generally two different ways to do this kind of integration (Press et al. 1996).

I. Averaged Monte Carlo integration
For function, f(x) in the [a, b] interval, by getting the averaged function, f , its integration can be obtained as follows: To calculate the averaged function, we first generate N random number with uniform distribution in the [a, b] interval and then we get the averaged function, f , according to Therefore, Eq. (71) can now be written as (see Fig. 1I)

II. Monte Carlo integrations, based on pair point method
In the first step, the random number,x i , is generated uniformly in the [a, b] interval [using Eq. (70)]. In the second step, the random number, y i , is generated in the [0, c] interval in which c ≥ f max and f max is the maximum value of the function f in the [a, b] interval. Therefore, we achieve pairs of numbers (x i , y i ) in a rectangle with b − a and c dimensions. By repeating the generation processes in the first and second steps N times, We should note that for negative values of y i , we decrease one unit from m and consider the absolute value of the function, f, to produce the related pair. This method is more appropriate for huge complex integrals and by increasing the number of repeating processes, N, we get a better solution. After specified iterations, we get to a converged solution and we would not need to add to the number of repetitions.

Laguerre function, xlag(N,y,nmax)
The function xlag(N,Y,nmax) will give us numerical values for the Laguerre function at each order. By accessing the Laguerre function at two successive orders n and n − 1, the Laguerre function at the order n + 1 will be obtained using the following relation (Arfken and Weber 2005):

Subroutine intp0(p0,ymin,ymax,ndat,nmax)
This subroutine will produce the Laguerre expansion coefficient of the splitting function, using Monte Carlo integration. The input of this subroutine is the splitting function and the output is an array which contains the subtraction of subsequent expansion coefficients and is denoted by p0.

Subroutine intR(R,ymin,ymax,ndat,nmax)
This subroutine will produce the expansion coefficients of a combined splitting function at the NLO approximation in each order, using the Monte Carlo integration. The input of this subroutine is the splitting function and the output is an array which contains the subtraction of the subsequent expansion coefficients and is denoted by rn.

Splitting functions, FPn0(n,x,nmax),…
The outputs of these functions are the numerical values of splitting functions which are multiplied by Laguerre coefficient in which the existing singularities are removed by the plus prescription method. The plus prescription takes advantage the following relation (Greiner et al. 1996):

Function E0Lag(y0,ELO,nmax)
Considering the Laguerre expansion for the evolution operator which was introduced in subsection "Laguerre expansion", the numerical values for this function can be obtained at each order of Laguerre expansion in terms of t and x variables.
Valance u quark densities in the LO approximation at energy scales Q 2 = 4, 50, 200 GeV 2 . Comparison with the CETQ4L and GRSV98LO parameterization groups has also been done Ghasempour Nesheli et al. SpringerPlus (2016) 5:1672 Functions qinq(z), … These functions are related to the parton densities at initial energy scale Q 2 0 = 2.56 GeV 2 and their combinations which could be found in Lai et al. (1997). These functions are calculated using the generated random numbers, based on the Monte Carlo program. The results of two CTEQ4L and CETQ4M fitting groups are used to give us the initial parton densities at LO and NLO approximations respectively.

Function Zeta(is)
This function is defined by Ellis et al. (1996)

Function S2(Y)
The S2(Y) function is defined by Ellis et al. (1996) where Li 2 (x) as a dilogarithm function can be approximated by: The other required functions and subroutines will be introduced in the respective sections.

The programs and the results
All programs are written in four parts:

Non-singlet sector at the LO approximation
In this case, we ate concerned with the evolution of valence quarks. So according to notations of section "DGLAP evolution equations", we just need to consider the DGLAP equation for q (−) i . To achieve this numerical solution, we do the following: 1. First, we should calculate numerically the expansion coefficients of splitting functions, using Eq. (37) where we choose the upper limit of summation in Eq. (36) equals to 30. To avoid the singularities which do exist in the splitting functions, we take into account the splitting function while we multiply the function with x rather than themselves and in the end, we display x times the parton densities. Therefore, the expansion coefficients are obtained via the following relation: in which (Furmanski and Petronzio 1980): x 3 3 2 + · · · for |x| ≤ 1. (80) The contribution of Dirac delta function in Eq. (81) is given by: Considering the rest of P (0) V (x), the final result for Eq. (80) would be: By applying the plus prescription,Eq. (76), the final result for P (0) n is given by: The integral in Eq. (84) is done numerically, using the pair point method [see Eq. (74)]. The result of integral is given by intp0 subroutine. The outputs of the program which contain the differences between two subsequent expansion coefficients will be saved in an array called p0(0:nmax) [see Eq. (39)]. 2. Now using Eq. (42) it is possible to calculate numerically the matrix A. The sum which is relating to A matrix in Eq. (42) is calculated in sumA(A,p0, nmax) subroutine, where p0 is the output of intp0 subroutine. Therefore A in sumA(A,p0,nmax) subroutine is a two dimensional array which is presented by A(0 : kmax , 0 : nmax) in our program.
Therefore, the sum in Eq. (42) can be calculated by doing three loops over I, k and n indices. The last index, I, contains by itself an additional sum. All these steps are performed in sumA(A, p0, nmax) subroutine. The general solution for the A matrix is as follows: 1. Now by considering the A matrix from step 2, and p0 from the first step, it is possible to numerically calculate the evolution operator at LO approximation, using Eq. (41).
The related program is given by ELOn(ELO,A,p0,nmax) subroutine in which the variable tat the LO approximation has been defined: . Ghasempour Nesheli et al. SpringerPlus (2016) 5:1672 The required quantities in the definition of t can be found in Gluck et al. (1998) [see Eq. (3)].
1. The convolution integral in Eq. (17) would now take the form: or alternatively: since in the CETQ4 parameterization (Lai et al. 1997), the initial densities are presented as x times the parton densities. Considering the Laguerre expansion of evolution operator, Eq. (89) can be written as: where E (0) n (t) is given by Eq. (41). How to use it in the numerical calculation was described in previous step. In Eq. (90) yq (−) i (y) represents the initial parton densities at energy scale Q 2 0 = 2.56 GeV 2 . The numerical solution of Eq. (90) is brought into the main part of program. For this propose we first need to perform a loop over x in the interval, say [0.001, 1] with the step 0.001 and then do the integration over y for each value of x. We perform the integration using the pair point method [Eq. (74)] which is specified in the program by the index pair. This program is run at three energy scales Q 2 = 4, 50, 200 GeV 2 where the wall clock time depends on the produced random numbers. The output of the program is labeled "Our result". The results are depicted in Figs. 2, 3 and compared with CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups.

Non-singlet sector at NLO approximation
In this case we should numerically solve Eq. (12) at the NLO approximation. We first need to calculate R − which was defined by Eq. (15). This section, like the previous one, can be divided to four parts.
1. The required splitting functions at NLO approximation are could be found in Furmanski and Petronzio (1980), Herrod and Wada (1980). Whenever it is needed we use the plus prescription technique to remove the singularities in the calculations [see Eq. (76)]. To obtain the Laguerre expansion coefficients, all the splitting functions should be multiplied by xL n (ln 1/x) and then we need to perform numerical integration as we did before [Eq. (80)]. It is required to separately perform the integration resulted from Dirac delta function in the splitting functions. These integrals together with the integrals whose singularities have been removed and also the rest of results, will produce the following functions which we need to run the program: xRnglag, xP0qqlag, xP1nglag, PF , PA, xPGlag, fPG, xPNFlag, fPNF . The subroutine which provides us with the Laguerre expansion coefficients of P − , R − is called intp0 (p0,rn,xmin,xmax,ndat,mmax). The contribution from Dirac delta function in the splitting function can be expressed by: The numerical values, resulted from Eqs. (91, 92) has been calculated in the intp0 subroutine. In continue, by adding all the contributions from the splitting function, the related Laguerre expansion coefficients can be calculated as: where P (1) − NS has been introduced in Ref. Furmanski and Petronzio (1980), Herrod and Wada (1980). The numerical solutions of Eq. (93) are calculated in the intp0 subroutine. The outputs of the subroutine are the differences between two subsequent coefficients of the expansion which will be put in p0(0:nmax) and rn(0:nmax) as two dimensional arrays. 2. Now by accessing the two arrays p0(0:nmax) and rn(0:nmax) it is possible to calculate the matrix A as before. The program does not need to be changed as we use sumA(A,p0,nmax) again and we simply change its inputs. 3. In this step, considering the matrix A which was made in the second step and the values of p0 and rn which were calculated in the first step, it is possible to calculate the expansion coefficients of evolution operators, E (0) n (t), E (1) n (t), E n (t) [see Eqs. (50, 51)]. The related program is presented in ENLOn(ENLO,A,p0,rn,nmax) subroutine, as it follows: At first the variable: should be used where the coupling constant at NLO approximation is given by:

The evolution operator at the NLO approximation in terms of the t variable is written as [see Eq. (50)]:
After performing all the numerical calculations in the mentioned subroutine, the output will be saved in a one dimensional array called ENLO.
1. This step is like step 4 of the previous part. The valence quark densities are obtained, using the following convolution integral: The expansion coefficients, E n (t), can be obtained via: Equation (99) is calculable with E0lag subroutine (see subsection "Function E0Lag(y0,ELO,nmax)") where E − is governed by Eq. (43) and since R − and P − are known, the expansion coefficients, E n (t), are calculable.
To perform the integration in Eq. (98), we follow the previous method. The only difference is in the initial parton densities which would be the ones mentioned in Lai et al. (1997)) but at NLO approximation. The results of running the programs for valence u and d quark densities are depicted in Figs. 4, 5 and compared with the results from the CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups.
To indicate the reliability of our numerical calculations, we should provide their statistical errors. This can be done, using the following relations: where "f " refers to related parton density. The last term in this equation indicated the error of calculations where 〈f 2 〉 and 〈f〉 2 are given respectively by: In these relations N denotes to number of points where we used in the Mote Carlo numerical integration.
As a result what we got for statistical errors of different valance densities at typical value Q 2 = 50 GeV 2 and at the NLO approximation are as following: xU v → %0.76 , xD v → %0.33. As can be seen these errors are <1 % which indicates enough reliability of our calculations for valence densities. The complete information for the other statistical errors at different energy scales for both LO and NLO approximations can be found in the appendix.

Singlet sector at the leading order
In the singlet sector, we encounter functions which possess a matrix form. According to subsection "DGLAP evolution equations", we should solve the Matrix evolution equation and also the equation for χ i , Eq. (14), in order to obtain the sea quark densities. First, the Gluon and q (+) densities are obtained. Then, using the q (+) and χ i , it is possible to get the q i (+) densities. And finally, having the valence and q (+) i , the sea quark densities for any separate flavor is extractable.
For this section, two separate programs have been written and we first take into account the solution for the gluon density.

Gluon distribution at leading order
To get the solution for gluon densities, we should note that in almost all parts of the calculations we would encounter matrix form. This section involves three steps.
1. First, we need to define the splitting function (Furmanski and Petronzio 1980;Herrod and Wada 1980) As can be seen, the splitting function is defined by a two dimensional array. As before the expansion coefficients of the splitting function should be determined which perform a three dimensional array: Equation (101) is in fact a matrix form which should be applied for all the elements in Eq. (100). It is obvious that we need to resort to the plus prescription technique to overcome the singularities of integration in Eq. (101). The concerned integration has been solved by the first method of numerical integration which is called the "averaged method" [see Eq. (73)]. The numerical solution of the related integral is inserted in the intp0 (p0,e1,e2,xmin,xmax,ndat,nmax) subroutine. The output of this subroutine is the difference between the two subsequent expansion coefficients which is put in a one dimensional array, named p0(0:nmax).
The projection operators e 1 and e 2 are finally given by Eq.  2. Secondly, we should computed the A and B matrices, given by Eq. (61). To define the related arrays, we first consider the upper index (k) and then the lower index (n) and then we proceed to the indices of the 2 × 2 projection matrices which are represented by IE, JE symbols. So A and B construct 4 × 4 dimensional arrays which can be represented by: The limit of the arrays are specified in the left hand side of Eq. (103). The IE, JE symbol is represented by 2 × 2 matrices. Equations (61, 63, 64) are three recurrence relationships. To get the related solutions, we first need to calculate them for n = 0. Therefore we will have: These are considered as the required constants in our calculations. To proceed, we would consider n = 1; therefore, we would have The above relations can be extended to n > 1. So we can acquire all the required values of A and B matrices. There are four A and B matrices which are given in the program by IE, JE(qq, qg, gq, gg). The general forms of the matrices are: (104) . Ghasempour Nesheli et al. SpringerPlus (2016) 5:1672 All the calculation are performed in the ABELO (E0,p0,e1,e2,nmax) subroutine. The inputs of this subroutine are e 1 , e 2 and p0 matrices and the outputs are the Laguerre expansion coefficients of the evolution operator, Eq. (59), at LO approximation, where they are put in the three dimensional array E0(n,IE,JE).
3. By achieving the Laguerre expansion coefficients, the evolution matrix operator in the Bjorken x space is calculable. Equation (55) is given by E0Lag function which which is called in the main program whenever is needed. We should note that at LO approximation, E ± are equal toE (0) . The combinations which give us the required densities are as follow: As before, the parton densities at initial energy scale Q 2 0 = 2.56 GeV 2 are taken form (Lai et al. 1997). The used density functions in the program are given by: The first function is related to q (+) . The second and the third ones are related to valence quarks. The rest are related to sea quarks except the last one which is for the gluon density.
The inputs of these functions are the generated random numbers by Monte Carlo method which produce numerical values for the related parton densities at initial energy scales. Since the density functions are appeared as xq, for q (+) we will have: Now, using the input density functions, the convolution integrals, Eq. (65) and the function E0Lag which is related to Eq. (55), it is possible to get the q (+) and gluon densities as in the following: (110) A, B(k, n) = 0 , for k > n.
qinq, qinU , qind, qinqUb, qindb, qinSb, qinG The integrals in Eqs. (117,118) are obtained numerically based on the "average method" [see Eq. (73)]. The above integrals are used in the main part of the program. The results include q (+) and gluon densities. We would require q (+) in further sections to obtain the sea quark densities. The results of running the programs for gluon densities at different energy scales are depicted in Fig. 6 compared to the results from the CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups.

Sea quark densities at LO approximation
This section contains two parts and for each part we require a separate program. At first, the equation which is related to χ distribution is solved. Then using the valence (obtained in 4.2) and q (+) densities which were obtained in subsection "Gluon distribution at leading order", we would be able to extract sea and gluon densities. This subsection is divided to two steps: 1. At this step we should first solve Eq. (9) for χ i distribution. The related solution is like the one for the non-singlet distribution (subsection "Non-singlet sector at NLO approximation") except that the replacement given by Eq. (53) should be carried out.
Therefore we will have [see Eq. (13)] where The analytical form for P (1)+ NS can be found in Petronzio (1980), Herrod andWada (1980). In the written program we need to define two functions (xRpolag, xP1polag) which should be replaced by previous ones in subsection "Non-singlet sector at NLO approximation" (xRnglag, xP1nglag). The other stages, are just the same as the ones in section "Non-singlet sector at NLO approximation". Since we have used R + and P + , the obtained coefficients are E + . And after that, it would be possible to obtain the χ i (x, t) distribution in Bjorken x-space by Eq. (18).
The initial χ i densities for separated quark flavor will be obtained, using the Eq. (66) as follows: (118) n,gq (t)L n ln y x yq (+) (y) + nmax n=0 E (0) n,gg (t)L n ln y x yG(y) dy y 2 . Comparison with CETQ4L and GRSV98LO parameterization groups has also been done Ghasempour Nesheli et al. SpringerPlus (2016) 5:1672 The above functions are added to the program by the following names: Having accesses to the initial parton densities from Lai et al. (1997), which are recognized by the tilde symbol, we will get the χ i densities for separate flavor as in the following: The integrals in Eq. (122) can be solved numerically, using the "average method" [see Eq. (73)]. The results would be the χ i distributions at different energy scales.

2.
It is now possible to get the sea quark densities at different energy scales, using Eq. (66): The required densities in Eqs. (123-125), including the valence (subsection "Non-singlet sector at NLO approximation"), χ i [Eq. (122)] and q (+) (subsection "Gluon distribution at leading order") distributions have been obtained before. The results of running the programs for sea quark densities at energy scales, Q 2 = 4, 50, 200 GeV 2 are depicted in Figs. 7, 8 and 9 and compared with the results from CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups.

Singlet sector at the next leading order
As before this section includes two subsections:

Gluon densities at NLO approximation
The required Laguerre expansion coefficients for the evolution operator is given by Eq. (97). The E (1) n (t) term in Eq. (97) is defined by Eqs. (68, 69). The only quantity in Eq. (69) which should be determined is R j . In the following, the gluon at NLO approximation can be calculated based on the following steps.
1. The calculation of R j is done by a subroutine called intR (R,xmin,xmax,ndat,nmax).
Evolution of gluon densities is done by Eq. (14) where P (1) (x) and R(x) in Eq. (16) are defined by: xkhid, xkhiS, qinq, qinU , qind, qinU , qinUb, qindb, qinSb, qinG . (122) Comparison with CETQ4L and GRSV98LO parameterization groups has also been done while matrix form for P (0) (x) is given by Eq. (100). As before, these matrices are defined by two dimensional arrays. The analytical expression for the splitting function at NLO approximation can be found in Furmanski and Petronzio (1980), Herrod and Wada (1980). The used functions in the program whose singularities have been removed by plus prescription technique [Eq. (76)] are: The required integrals are done numerically in the intR subroutine. The output of the program is put in three dimensional arrays which are called R(0:nmax,2,2). 2. At this step the evolution operator E (0) should be calculated; this has been done in subsection "Gluon distribution at leading order". The only difference is in the definition of the t parameter which should be redefined at NLO approximation. The matrices A and B which are used to get the evolution operator are like before. So in this step we can use a program similar to what has been written in subsection "Gluon distribution at leading order". This program contains the following functions and subroutine: After computing A and B matrices, as before we will have In this equation, we use different values for t and we denote Eq. (128) by the function: E0t (nP,IP,JP,t,A,B,nmax). The inputs of this function are the indices of the 2 × 2 matrix elements, A and B matrices and the t variable. The output of the function is the numerical values of Laguerre expansion coefficients for evolution operator. 3. To get Ẽ (1) n in Eq. (68) we should solve the integral in Eq. (69). For this purpose we should first solve the sum in the integrand of this equation. So Eq. (69) can be written as: where The E (0) n (t) terms can be calculated using the E0t function [see Eq. (128)]. The R j was calculated before as the R(0:nmax,2,2) array (step.1) which can be considered as the input of the program. So we should write a program which calculates the sum in Eq. (130) by taking into account the Dirac delta function and also do the matrix multiplication. The details of the program can be found in the following subroutine: The outputs are the four elements of the matrix which are related to the sum in Eq. (130). 4. In this step we first calculate the integral in Eq. (129) and we would obtain the related Laguerre expansion. Then we do the loop over the order of Laguerre expansion since we are going to calculate the integral in Eq. (129) at each order n. The τ variable is obtained by generating the random number, RAN3, in the interval [0 , t NLO ] so as Press et al. (1996) Now we give the random number to the SUMERE subroutine. Then we call this subroutine for each generated random number in which we are able to calculate the sum in Eq. (130). To increase the precision of calculations we need to generate more random numbers which consequently increase the wall clock time. The integrals can be calculated, using the "average method" [see Eq. (73)]. We should note that four integrals are calculated which are in fact the elements of the matrix evolution operator Ẽ (1) n (t). The solutions of these integrals are put in the three dimensional array Et1(−2:nmax,2,2) where the first two terms of this array for all elements of the 2 × 2 matrices are zero [see Eq. (68)]. We continue to calculate the required expression in Eq. (68). The E (1) n (t) term is appeared as an array by the name E1(0:nmax,2,2). This array is put inside the n th loop. The term, relating to n = 0 is computed, using the first two terms which were introduced in Eq. (68). The other terms, relating to the Laguerre expansion at NLO approximation are obtained by iterating the loop over n, IE, JE indices. In the end, using this subroutine and the function E0t which was obtained before, we can calculate the related Laguerre expansion coefficients at the NLO approximation, based on Eq. (97). All the required mentioned tasks are gathered in a subroutine called intE (En,A,B,R,ndat,nmax). This is the most important part of the program. The outputs are the expansion coefficients of the evolution operator as 2 × 2 matrices. 5. Now by getting the Laguerre expansion for the evolution operators and the parton densities at the initial energy scale Q 0 2 = 2.56 GeV 2 , we can obtain the evolved parton densities at any desired energy scale. The process is like the one for the LO approximation. In this case we will have The initial parton densities at the NLO approximations are taken from Lai et al. (1997) which are designated in Eqs. (132, 133) by tilde symbol. The list of the required functions is: The integrals in Eqs. (132, 133) are numerically calculated, using the "average method" [see Eq. (73)]. The solution of integrals is brought into the main part of the program. The results are gluons and q + densities (singlet sector); the singlet densities will be used in the next section to obtain the sea densities. The outputs of this program are the evolved gluon densities at different energy scales which are depicted in Fig. 10 and qinq, qinU , qind, qinqUb, qindb, qinSb, qinG are compared with the CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups. If we wish to calculate the statistical error for gluon density, we should resort to Eq. (99-a). What we get for the required error at the typical energy scale Q 2 = 50 GeV 2 in the NLO approximation is %1.78 which indicates good precision in our calculations for gluon densities (for more information, see the "Appendix").

Sea quark densities at the NLO approximation
The objective of this subsection is to obtain the sea quark densities at the NLO approximation. The entire procedure is like what we have been done for LO approximation (see subsection "Sea quark densities at LO approximation"). The required relations to obtain the sea quark densities are given by Eqs. (123-125). Here, two programs should be run. The first program gives the gluons [see Eqs. (132,133)]. The second program, which uses the results of the first program (q (+) ) and the results of the program which gave us the χ i densities [see Eq. (122) (123)(124)(125)] at different energy scales in the NLO approximation. The results and the comparisons with the CTEQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups are presented in Figs. 11, 12 and 13 for different quark flavors.
As can be seen from Fig. 13 our results for the strange sea quark densities are in good agreement with the CETQ4 M and GRSV98NLO parameterization groups. The results from Ref. Coriano and Savkli (1999) indicate completely different behavior with respect to the fitting parameterization models as well as with respect to our results. This confirms the validity of our calculations for numerically obtaining the evolved parton densities.
To provide the statistical error for the sea parton densities we need again to resort to Eq. (99-a). The obtained errors at the typical energy scale Q 2 = 50 GeV 2 for the NLO approximation are as following (see appendix as well): Once again the small values for the statistical errors indicated enough precision of the employed numerical integration to evolve the parton densities

Conclusion
In this paper, we have presented numerical solutions for the DGLAP evolution equations, based on the Laguerre polynomials expansion (Furmanski and Petroznio 1982a, b). Although people can use other methods especially in the Mellin moment space, the method which we used in this article has this specific feature that we do not need to change the space of calculations. In fact, all the computations have been done in the Bjorken x-space. We have tried to explain all the steps of performing the FORTRAN codes which produce parton densities at high energy scales. Since we have just used FORTRAN package, it means that all calculation have been done numerically. The main program can be requested from the authors via the E-mail address: a.mirjalili@yazd.ac.ir. The results are in good agreement with CETQ (Lai et al. 1997) and GRV (Gluck et al. 1995) parameterization groups. This confirms the validity of our numerical solutions for the DGLAP evolution equations. Our results for parton densities are much better xŪ → %0.17 , xD → %0.2 , xS → %0.076. Fig. 10 Gluon densities in the NLO approximation at energy scales Q 2 = 4, 50, 200 GeV 2 . Comparison with CETQ4M and GRSV98NLO parameterization groups and the results from Ref. Coriano and Savkli (1999) has also been done. Our results for gluon densities are in better agreement with the CETQ4M and GRSV98NLO results rather than the results from Ref. Coriano and Savkli (1999) Fig. 11 Sea u quark densities in the NLO approximation at energy scales Q 2 = 4, 50, 200 GeV 2 . Comparison with the CETQ4M and GRSV98NLO parameterization groups has also been done Ghasempour Nesheli et al. SpringerPlus (2016)