Analytical solutions of the QCD$\otimes$QED DGLAP evolution equations based on the Mellin transform technique in NLO approximation

In this paper we present a new and efficient analytical solutions for evolving the QED$\otimes$QCD DGLAP evolution equations in mellin space and obtain the parton distribution functions (PDFs) in perturbative QCD including the QED corrections. The validity of our analytical solutions, which have done in the next to leading order QCD and the leading order QED approximations, are checked with the initial parton distributions from newly released CT14QED global analysis code (Phys. Rev.D93,114015(2016)). The evolved parton distributions are in good agreement with CT14QED code and also with those from APFEL (Computer Physics Communications 185, 1647 (2014)) program.


Introduction
The new search at the LHC demands knowledge of the photon distribution function at different values of x and Q 2 . The well-known Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) integrate-differential evolution equations give the parton distribution functions in the pertobative QCD [1,2,3,4]. There exists many published papers attempting to give solutions of DGLAP equations in the realm of QCD, mostly 5 based on the global parameterization of PDF's, but only a few have also included corrections coming from QED. The MRST group [5,6] being the first known group to that effect.
Recently the NNPDF collaboration [7,8] and CT14QED group [9] reported the new results on this issue. A precise knowledge of the parton distribution functions of the proton is presented in Refs. [10,11] in order to make predictions for the Standard Model and beyond the Standard Model processes at hadron

Evolution Equations
The singlet parton distribution functions f i (x, Q 2 ) obey the DGLAP evolution equations with QED corrections [25] in x space, as P 11 P 12 P 13 P 14 P 21 P 22 P 23 P 24 P 31 P 32 P 33 P 34 and the DGLAP evolution equations with QED corrections for the non-singlet parton distribution functions are as follows, ∂f i ∂logQ 2 = P ii ⊗ f i i = 5, . . . , 9 (2) where P ij and P ii are the splitting functions and are given in Ref. [20] with details, and ⊗ denotes the convolution integral We utilize a PDF basis for the QCD⊗QED DGLAP evolution equations. This basis define by the following singlet and non-singlet PDF combinations [26], u +ū + c +c + d +d + s +s + b +b Here in the next subsections we bring out the solutions of the QCD⊗QED DGLAP evolution equations with more details. Our solutions is done in the next to leading order QCD and the leading order QED approximations. 40 In the singlet sector, f in equation (4) is a four component vector. The Mellin transform is defined as,

The singlet PDFs with LO QED corrections
Therefore, the evolution equation (1) factorizes into the form where P N (Q 2 ) is the 4 × 4 splitting matrix. We can expand this splitting matrix in the next-to-leading order QCD and the leading order QED approximations as follows, P N (Q 2 ) = a s (Q 2 )P (1,0) 0 (N ) + a 2 s (Q 2 )P (2,0) (N ) + a(Q 2 )P (0,1) 0 where a s ≡ αs 4π and a ≡ α 4π are the QCD and QED running couplings, respectively. We can separated the splitting matrix into two parts of QCD and QED. Then we have, The QED splitting matrix, P (0,1) 0 (N ), and the QCD splitting matrices P (1,0) 0 (N ) and P (2,0) (N ) are rep-45 resented as follows, The n u and n d parameters are the number of up and down-type active quark flavors, respectively, and where, n f is the number of active flavors. The moments of Altraelli-Parisi function may be found in Refs. [26,27,28,29,30,31,32,33,34,35], noting that P z i = +γ (i)z /2. The running couplings 50 are given by where, β (αs) 0 where, e i is the relative electric charge of the quarks and N c denotes the multiplicity due to color degrees of freedom, i.e.,N c = 3 for quarks. The n l parameter is the number of active lepton flavour.
The evolution equations of these couplings are given by the QCD and QED β-functions as where, the beta functions for the strong and electromagnetic couplings are expanded to the appropriate order, as We obtain the solution of the equation (7) in a compact form, where E N (Q 2 , Q 2 0 ) is the evolution matrix, as where E N QCD (Q 2 , Q 2 0 ) and E N QED (Q 2 , Q 2 0 ), are the QCD and QED evolution matrices, respectively. The QCD evolution matrix, E N QCD (Q 2 , Q 2 0 ), is the solution of following equation Defining t as, The equation. (17) can be write as follows, where, R = P (2,0) (N ) − Then the full solution of equation (19) can be write as a power series of α s as, Substituting the equation (21) into equation (19), we obtain the following equation for the new evolution matrix U : where, U, P  Using the projection operators, we can write-down P where, the projection operators are given by We have where, 1 is the unit 4 × 4 matrix. The solution of the equation (20) can now be written as Since e 1 + e 2 + e 3 + e 4 = 1, we have an obvious identity matrix as, U = e 1 U e 1 + e 1 U e 2 + e 1 U e 3 + e 1 U e 4 +e 2 U e 1 + e 2 U e 2 + e 2 U e 3 + e 2 U e 4 +e 3 U e 1 + e 3 U e 2 + e 3 U e 3 + e 3 U e 4 Inserting the Eq. (23) into Eq. (22), yields Therefore, The final solution for the QCD evolution matrix, E N QCD (t), in the next to leading order approximation can be obtained from the equation (21) as where, a s0 is a s at a initial scale of Q 2 0 . Now, the leading order QED evolution matrix, E N QED (a, a 0 ) in terms of a and a 0 at the scales of Q 2 and Q 2 0 , respectively, can be represented as The matricesẽ 1 ,ẽ 2 ,ẽ 3 andẽ 4 denote the corresponding projectors, where,λ 1 ,λ 2 ,λ 3 andλ 4 are the eigenvalues of QED matrix, P With the well-known inverse Mellin transform [36] the parton distribution functions can be derived in x space, where the contour of the integration lies on the right of all singularities of M k (N = c + ωe iφ , Q 2 ) in the complex N-plane. 80

The non-singlet PDFs with LO QED corrections
For the non-singlet distribution functions, all of the parameters in Eq. (2) are scalar. Since the Mellin transformation turns convolutions into products, the evolution equation for QCD part, in the next to leading order approximation, becomes With a change of variable a s instead of logQ 2 in equation (31), we have Now we expand the right-hand side of equation (32) to get Eq. (32) and Eq. (33) are formally solved by introducing an evolution operator, E N QCD (Q 2 , Q 2 0 ). In the non-singlet case this evolution operator is just a scalar function of N. This evolution operator satisfies the equation (33), then we have Solution of this evolution equation is straightforward in the non-singlet case; expanding beyond the leading terms, we obtain [37,38,39,40,41] With the same procedure we have the QED evolution operator as, Finally, we simply have the non-singlet PDFs including QED corrections as The equation (38) has a boundary condition, E N (Q 2 0 , Q 2 0 ) = 1, that renders it an entirely perturbation object. Now by using this method we obtain operator E N (Q 2 , Q 2 0 ) for all of the distribution functions in equation (2). Our results for these functions presented in Table. (1).

The impact of NLO QED corrections 85
In this section, the discussion starts by generalizing the present analytical method for the singlet and non-singlet evolution equations with NLO QED corrections.

The singlet PDFs with NLO QED corrections
Here, we study the impact of the NLO QED corrections to the singlet DGLAP evolution equations.
These equations with NLO QED and QCD corrections include mixed terms of α s and α. Therefore, the QCD and QED β-function are proportional to α s and α. The evolution equations of couplings a s = αs 4π and a = α 4π are given by where the beta functions can be expanded to the appropriate order, as follows where Now, the splitting function matrix that determines the evolution of the singlet PDFs in NLO QCD and NLO QED approximations are as follows where the pure QCD term is given bỹ and the terms including the QED coupling are given bȳ The evolution kernels P (1,1) and P (0,2) are terms of order O(a s a) and O(a 2 ), respectively, that add to the DGLAP evolution equations. These evolution matrices are given by The splitting functions corresponding to these evolution matrices can be found in Refs. [12,13].
In the following, for calculating the NLO QCD evolution matrix, E N QCD (a s , a s0 ), we used the method represented in subsection (2.1). Now, we consider the series expansion for splitting functions Eq. (41) and for the beta function Eq. (39) and sort everything in a power series in a s and a, yields Notice that in order to solve the above equation, we expand a in terms of a s . Then we can choose a = A + Ba s + O(a 2 s ). Therefore, we can rewrite Eq. (46) as follows Introducing futher the abbreviations Finally, the solution of equation (47) at NLO QCD and QED approximations with the method represent in subsection(2.1) is given by The QED part of evolution matrix satisfy in the following equation, where, thet parameter is 1 a(Q 2 ) . In this case, we expand a s in terms of a. Therefore we choose a s = Ca+O(a 2 ). We ignore the last terms in this expansion, because we work in NLO QED approximation.
We use this expansion and substitute it in equation (50), then we have Now, in the same method with QCD part, we define the new following parameter Therefore, the QED evolution matrix obtain in the following form, We calculated the values of fixed parameters A, B and C for the different active flavour numbers. We listed them in Table (2).

The non-singlet PDFs with NLO QED corrections
In the non-singlet case, the evolution equations are scalar, so that the commutation relations in Eq. (22) are vanished. In the same procedure established in subsection (2.2), we calculate the evolution operators for QCD part in the next to leading order QCD and QED approximations for the non-singlet case.

95
These QCD evolution operators for the non-singlet PDFs introduced in Eq. (5) are listed in Table (3).
With the same procedure, we can obtain the QED evolution operators for the non-singlet PDFs in 100 Eq. (5) with NLO QED corrections. We give these operators in table (4).

Test the accuracy of the solutions
In the section (2), we analytically solved the QCD⊗QED DGLAP evolution equations in NLO QCD and LO QED approximations. We obtain the singlet and non-singlet parton distribution functions with QED corrections in x space and for Q 2 > Q 2 0 (GeV 2 ).

105
The main question here is: How we can check our analytical solutions? The best way to be sure about the correctness of these solutions is to take the initial PDFs from a public code such as the CT14QED global analysis code [9] at a fixed scale of Q 2 0 and then we can calculate all of the parton distribution functions inside proton at some values of Q 2 . In order to do this, we run the CT14QED program [9] with P γ 0 ≤ 14% for the inelastic photon PDF in initial scale of Q 0 = 1.295GeV . Then, we evolved this initial Comparison between the APFEL (NNPDF2.3QED) predictions and the CT14QED evolution with our 120 results for the valance quarks at different values of Q 2 is shown in Fig. (1). An excellent agreement is found for all flavors. It is clear from Fig. (1) that increasing the value of Q 2 would lead to a decrease in the value of decline rate the quark valence distribution functions. This means that the presence of photon distribution function influences on the decline rate.
The sea quarks, gluon and photon distribution functions are plotted in Fig. (2), where PDFs have 125 been evolved from Q 0 = 1.295GeV up to Q = 10 3 GeV and they are compared with the CT14QED and APFEL (NNPDF2.3QED) PDF sets. It is obvious that a good agreement is achieved.
One observation from Fig. (2) is that the sea-quark distributions , especially c and b quarks distributions, have more impact on the photon PDF than does the initial photon distribution at high Q 2 .
These photon distribution become more significant at high Q 2 where more photon are produced through 130 radiation of the quarks. It is worth to notice that, the photon distribution functions are larger than the sea quarks distribution functions at high scale of energy and for the large values of x.
In tables (2) and (3) we present the Patron Evolution program [26] results and the CT14QED code In conclusion, we found a good level of agreement for all of the comparison performed in this section.
These guarantees that we implements correctly the QCD⊗QED evolutions, therefore it can be used in  any global parametrization of PDFs with experimental data.
Our analytical solutions can be used to interpret the LHC data in a global fit parameterization.
To do this, recent ATLAS measurement of high-mass Drell-yan dilepton production data combined with inclusive deep-inelastic scattering cross-section data can perform a solid determination of the proton PDFs including the photon PDF [14]. 150

Summary and Conclusions
The analytical solution is performed based on the Mellin transforms to obtain the QCD⊗QED parton distribution functions. Our calculations are done in the NLO QCD and NLO QED approximations.
To be sure about the correctness of analytical solutions, we take the initial PDFs from the newly release CT14QED global parameterization at a fixed scale of Q 0 = 1.295GeV . Therefore, here we only checked   3. In our solutions, the QED running coupling is not constant and Q 2 dependence of this parameter is considered.
Last, but not least, these analytical solutions can be used to determine the QCD⊗QED parton distribution 170 functions via a global parameterization to the experimental data in the future.