Analytic derivation of the photon structure function $$F^{\gamma }_{2} (x,Q^{2})$$ at the leading order and the next-to-leading order approximations based on the Laplace transform method

In this paper, we present the analytical solutions, based on the Laplace transform method, for the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi evolution equations of the photon at the leading order (LO) and next-to-leading order (NLO) approximations in perturbative QCD. Using these solutions, we derive the singlet, non-singlet, gluon, and photon distribution functions of the photon and also the photon structure function $$F^{\gamma }_{2} (x,Q^{2})$$F2γ(x,Q2) at the LO and NLO approximations. We show that the resulting distribution functions are in agreement with the results of the parameterization model formulated by Gluck, Reya, and Vogt (GRV) (Phys Rev D 46(5):1973, 1992). Moreover, our numerical results of $$F^{\gamma }_{2} (x,Q^{2})$$F2γ(x,Q2) are comparable with the results achieved by Aurenche, Fontannaz, and Guillet (AFG) (Eur Phys J C Part Fields 44(3):395–409, 2005) and also with the experimental data released by the L3, DELPHI, OPAL, ALEPH and PLUTO collaborations.


Introduction
In the electromagnetic (EM) interactions, the photon mediates the EM force between EM charged objects. If the photon, in these kinds of interactions, fluctuates into a charged fermion-antifermion pair and one of them interacts with a gauge boson, then the parton structure of the photon is revealed. When fermions are a lepton-antilepton pair, the interaction process can be calculated in the framework QED, but if they are a quark-antiquark pair, this process is described in the framework of QCD. The calculations in the framework of QCD are complicated, because the spectrum of fluctuations is richer. In the photoproduction process, therefore, both the QCD and the QED contributions have to be taken into account in the calculations. There are two types of a e-mail: zarrin@phys.usb.ac.ir (corresponding author) photons in the interactions: the direct photon and the resolved photon. The direct photon takes part in hard interactions and its structure is not revealed, but the resolved photon first fluctuates into a hadronic state and its structure is revealed. The resolved photon includes two contributions, point-like and hadronic-like, which must be taken into account in the calculations of the structure function and the cross section [1]. In deep inelastic scattering, the results of the photon structure function F γ 2 (x, Q 2 ) have been reported by the large electronpositron collider (LEP) [2][3][4][5][6][7] and the electron-proton collider HERA [8][9][10].
After observing the photon structure function in different interactions, theorists have been interested in studying such a structure function. Since the photoproduction processes are susceptible to the properties of the photon itself, the electron-positron and the electron-proton colliders at low energy regions can be used as useful tools for probing the photon structure like that of the proton. On the other hand, there are many differences between the hadronic structure function and the photon structure function. For example, at low Q 2 , the hadronic structure function decreases but the photon structure function grows even without considering the QCD corrections [1,11].
So far, the photon structure function has been suggested by various theoretical methods: Aurenche et al. [12] obtained the photon structure function beyond the leading logarithm approximation. They modified the hadronic component by taking into account the MS factorization scheme. By fitting all the available data on the photon structure function, Gordon and Storrow [13] determined the parton distribution functions of the photon. In fact, they attempted to update these functions by constraining the input gluon distributions through fitting them to TRISTAN jet data; Vieira and Storrow [14] obtained the quark distribution functions of the photon and the photon structure function using the ansatz of Rossi [15] generating regularized asymptotic solutions to the inhomoge-neous Altarelli-Parisi equations; Kolanoski and Zerwas [16] calculated the total cross section for hadron production by two photons. They presented the photon structure function at LO and NLO by the quark parton model and the vector meson dominance model; Berger [17] calculated the photon structure function at LO and NLO based on the operator product expansion and gave a unified picture of the hadronic and point-like pieces. He obtained the point-like contribution of the photon structure function in the MS scheme and the hadronic-like contribution from the vector meson dominance model.
To obtain the parton distribution functions at higher energies than the initial scale Q 2 0 , one can use the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [18,19]. The DGLAP evolution equations are a set of coupled integro-differential equations that give the parton distribution functions in perturbative QCD.
Recently, the QCD DGLAP evolution equations of the proton at LO and NLO have been solved by using the Laplace transform [20][21][22][23] and the Mellin transform methods [24,25] and have consequently yielded the proton structure function and singlet, non-singlet, and gluon distribution functions inside the proton. Furthermore, utilizing these methods, one can solve the QCD⊗QED DGLAP equations at LO and NLO [26,27]. In these equations, the QED corrections are considered in the QCD calculations. Such corrections are critical and affect the distribution and structure functions and consequently modify the cross section of processes. In Refs. [26,27], it has been shown that the photon contribution at high x is comparable with t and b quark contributions in the proton [26,28].
In this paper, we apply the Laplace transform method to solve the QCD⊗QED evolution equations of the photon at LO and NLO QCD and also at LO QED. On this basis, we extract the singlet, gluon, and photon distribution functions of the photon at the scale Q 2 as follows: where , and M γ 0 (x) are the singlet, non-singlet, gluon, and photon distribution functions at the initial scale Q 2 0 , respectively. In the above equations, the functions F, G and M can be obtained by using the splitting functions P qq (x, Q 2 ), P qg (x, Q 2 ), P gg (x, Q 2 ), P qγ (x, Q 2 ), P γ q (x, Q 2 ), and P γ γ (x, Q 2 ) and the distribution functions at the initial scale. Moreover, using the singlet and non-singlet distribution functions, one can get the photon structure function at different scales Q 2 .
The rest of the present paper is organized as follows: In Sect. 2, we present general solutions for decoupling the QCD⊗QED DGLAP evolution equations of the photon at LO in both QCD and QED by applying the Laplace transform method. In Sect. 3, we use the same method to calculate the QCD⊗QED DGLAP evolution equations at NLO QCD and LO QED. In Sect. 4, we also present an analytical solution of the non-singlet sector of the photon structure function based on the Laplace transform method. In Sect. 5, the numerical results of the photon structure function, obtained by the Laplace transform method, are compared with the available L3 [2,3], OPAL [4,5], DELPHI [6], ALEPH [7] and PLUTO [29] data as well as the GRV [30] and AFG [31] results. Also in this section, we compare the results of the singlet, non-singlet, and gluon distribution functions of the photon with the GRV results [30]. At the end of this section, we present the photon distribution function of the photon and show a comparison of this distribution function with the s-quark and c-quark distribution functions at NLO adopted from the GRV results. Appendix A gives a brief explanation of the Laplace and the Mellin transforms. In Appendix B, we present the functions of a set of coefficients in the Laplace s space, which are dependent on the splitting functions.

Master formula for the LO corrections
The QCD⊗QED DGLAP evolution equations of the photon for the parton density functions of the quark, antiquark, gluon, and photon obey the following inhomogeneous equations, respectively [1,32]: These equations include the QCD corrections in the QED calculations. In the above equations, q γ i ,q γ i , g γ , and γ γ represent the quark, antiquark, gluon, and photon parton density functions of the photon, respectively. The P i, j (x) are the Altarelli-Parisi splitting kernels at one loop correction, which are defined as follows: where n f denotes the number of active massless flavors, C F = 4 3 , and T R = 1 2 . It should be noted that, at LO, P gγ = P γ g = 0 due to the missing photon-gluon coupling, but these splitting functions at NLO are not zero [33]. By summing Eqs. (4) and (5) and using F where q k − e 2 x (q k +q k ) are the singlet, gluon, non-singlet, and photon distribution functions, respectively. Here, A = n f k=1 e 2 q k and e 2 = 1 n f n f k=1 e 2 q k . In Eqs. (9)-(11), the symbol ⊗ represents the convolution integral given in Appendix A.
Equations (9)-(11) have the same structure as those mentioned in Ref. [26]. However, in the present paper, there are three main differences between our calculations and those done in Ref. [26]: First, in our work, the contributions ofP qq andPqq are zero in the evolution functions of the photon. Second, the non-singlet contribution of the photon structure function is considered here in solving these equations, while it has not been included in Ref. [26]. Third, in Ref. [26], the proton structure function has been obtained only at high Q 2 and x smaller than 0.1, whereas, here, the photon structure function is computable for all of x and Q 2 ranges.
As shown in Refs. [20][21][22], the Laplace transform L of the convolution factors is simply the ordinary product of the Laplace space of the factors. Thus, in order to make Eqs. (9)-(11) easily solvable, we transform these equations to the Laplace space. To do this, we apply the variable changes x≡ exp(−v) and y≡ exp(−ω) together with the Laplace transform (a brief explanation of this transform and its reverse is given in Appendix A), as a result of which we have where we used the following notations: and then the Laplace transform converts Eqs. (9)-(11) into three coupled ordinary first-order differential equations in s space which can be written as where the coefficients L O , L O , L O , and ϒ L O are the LO splitting functions in the Laplace space s [22]: where ψ is the digamma function and γ E is the Euler constant. By introducing the new variable τ as ∂τ (Q 2 ) ∂ ln(Q 2 ) = α s (Q 2 ) 4π , the coupled first-order differential Eqs. (14)- (16) in s space can be rewritten as To make calculations easier at the LO and NLO approximations, we consider the following expressions [20,26]: where the constants a i j and b i j can be obtained by the fitting functions [26]. The percentage relative errors of Eqs. (26)-(28) are given in Table 1. As can be seen, the maximum value of these percentages is less than 2.5%. The coefficients a i j and b i j in Eqs. (26)-(28) are determined within a few percent accuracy.
To solve the QCD⊗QED evolution equations at LO, we need to transform Eqs. (23)-(25) a second time, this time from τ space to u space. Therefore, in u space and s space, we have where the Laplace transforms are as follows: To solve Eqs. (29)-(31), we consider the parameter a 11 as an expansion parameter. On this basis, the solutions of Eqs. (29)-(31) can be obtained as three power series of the parameter a 11 . Since this parameter is much smaller than one, we use these power series only up to the second sentence. For example, in the interval 3 ≤ Q 2 < 20, we have a 11 = 0.0109. Then the second sentence is about two orders of magnitude larger than the third sentence. By setting a 11 = 0 in Eqs. (29)-(31), the first approximation of the function α e (Q 2 ) One can easily solve Eqs. (33)- (35) and obtain the functions F 1 , G 1 , and M 1 as follows: Then, by repeating the above processes, we have where F 2 , G 2 , and M 2 are the distribution functions of the photon in s and u spaces at LO in the second approximation. Also, f 0 and m 0 are Similar to the first approximation (a 11 = 0), the functions F 2 , G 2 and M 2 are easily obtained as follows: where A i , B i , and C i (with i = 1...5) are given in Appendix B.

Master formula for the NLO corrections
In this section, we intend to extend our calculations to the NLO approximation for the distribution functions of the photon. The QCD⊗QED evolution equations of the photon at NLO QCD can be written as Eur. Phys. J. C (2020) 80:319 To solve Eqs. (47)-(49), again we need the Laplace trans- to the s and u spaces, respectively. In this regard, the solutions of the evolution equations of the photon at NLO can be converted to  [34]. Since a 21 and a 31 are smaller than one, they can be used as expansion parameters. Therefore, by setting a 21 = 0 and a 31 = 0 in Eqs. (50)-(52), these equations can be rewritten at NLO as follows: By solving Eqs. (53)-(55), one can obtain the gluon G 1 , singlet F 1 , and photon M 1 distribution functions in the first approximation at NLO:  58) to v and τ spaces. In fact, they are just used in the second approximation as an initial function. To obtain the next approximation (a 21 = 0, a 31 = 0) of the functions F, G, and M, we use the transforms F 1 → F 2 , G 1 → G 2 , and M 1 → M 2 and also replace f 0 , g 0 , and m 0 with f 0 , g 0 , and m 0 in Eqs. (53)-(55), respectively. Thus, we have where Eur. Phys. J. C (2020) 80:319 To get the final equations, we have to return them to the usual spaces, namely the v and τ spaces. For this purpose, we use the Laplace inverse transform and the convolution integral (see Appendix A). Therefore, the distribution functions are obtained as follows: where the kernels K i j (with i = f, g, γ and j = 1, 2, 3, 4) in τ and v spaces are To get the numerical results of Eqs. (68)-(70), we need to know the distribution functions at the initial scale. For this aim, we use the fitted results of GRV [30] at the initial scale , in a good approximation, we take the photon starting distributions as in the following equation: Note that the above equation for the photon distribution function of the photon is similar to the photon distribution function of the proton which is obtained in Ref. [28].

The non-singlet sector
In this section, we provide the solution of the non-singlet distribution function by using the Laplace transform method at LO QED and LO QCD as well as at LO QED and NLO QCD. The QED⊗QCD evolution equation for the non-singlet density function of the photon is as follows [32]: where the kernel P qqns is the non-singlet splitting function and a(x) describes the γ → qq splitting. By using F γ ns = n f k=1 e 2 q k − e 2 x [q k +q k ], as introduced in Ref. [1] and Eq. (73), one can write the non-singlet distribution functions at LO and NLO as follows: 4π where the kernel Z ns (x) is Using the convolution integral, the relation ∂τ ( and also the change variables x ≡ exp(−v) and y ≡ exp(−w), we can convert Eqs. (74) and (75) as follows: To solve Eqs. (77) and (78), we need to use the Laplace transform from the v and τ spaces to the s and u spaces, respectively. On this basis, Eq. (77) is converted to a firstorder differential equation, whose solution is as follows: Similarly, Eq. (78) is also converted to an ordinary firstorder differential equation, whose solution in the first approx- and in the second approximation Note that to obtain Eqs. (80) and (81), we use Eq. (28). Taking all the above into consideration, the non-singlet distribution function F γ ns (v, τ ) at LO and NLO can be obtained: where the kernel K i jns = L −1 L −1 k i jns (s, u); τ ; v and the kernels k i jns (s, τ ) are presented in Appendix B.
Eur. Phys. J. C (2020) 80:319 For the numerical investigation of Eq. (82), we apply the non-singlet distribution of the photon at the initial scale, F γ ns (x, Q 2 0 ), which is obtained from the fitted results of GRV [30]. Finally, recalling that v ≡ ln(1/x) and τ = 1 4π and Eq. (82) back into the usual spaces, i.e. the Bjorken-x and virtuality Q 2 .

Results and conclusion
In the previous sections, we obtained the singlet, non-singlet, gluon, and photon distribution functions of the photon at the LO and NLO approximations by using the Laplace transform method. In this section, we shall present the numerical results which can be extracted from these functions. For the verification of our solutions, we compare our numerical results of the photon structure function with L3, OPAL, DELPHI, PLUTO, and ALEPH data and also our numerical results of the different kinds of distribution functions of the photon at the arbitrary scale Q 2 with those from the GRV results.
The results of Eq. (68) are depicted in Fig. 1. It indicates the evolution of the singlet distribution function of the photon at scales Q 2 = 10.8, 12, 15.3, 76.4, 120, and 780 GeV 2 . In this figure, the solid and dash lines present the solution which is resulted from the Laplace transform method at NLO and LO, respectively, and also the dot lines show the GRV results [30] at LO and NLO. In order for the difference between our results and the GRV results to be more visible, we display the ratio F γ s (G RV ) F γ s (T heor y) for different scales at LO and NLO in Fig. 2. As can be seen in this figure, the maximum ratio at high x is about 20%. Note that all the figures discussed in this section, the colored bands represent the error caused by Eqs. (26)- (28). Figure 3 shows the results of the gluon distribution function G(x, Q 2 ) = xg(x, Q 2 ) of the photon and their comparison with the LO and NLO analyses of the GRV results [30] at scales Q 2 = 10, 20, 100, and 1000 GeV 2 . Figure 4 displays the ratio G γ (G RV ) G γ (T heor y) for different scales. As a numerical illustration for our analytical approach to the study of the photon structure function F γ 2 (x, Q 2 ) at the LO and NLO approximations, we compared our results with the photon structure function obtained from GRV and depicted them in Fig. 7. Also, a comparison with the LEP-L3 data [2,3] at Q 2 = 10.8, 15.3, and 120 GeV 2 has been shown there. In the LEP-L3 data, the total systematic error is calculated from the quadratic sum of the systematic error and the additional systematic error due to the dependence on the Monte Carlo model. Moreover, in Fig. 7, we showed a comparison of the photon structure function with the LEP-OPAL data [4,5] at scales Q 2 = 10.8, 12, 76.4, and 780 GeV 2 . Note that the OPAL data is only calculated with the Eur. Phys. J. C (2020) 80:319 Fig. 8 A comparison of the photon structure function F γ 2 at different initial scales Q 2 0 = 1GeV 2 and Q 2 0 = 0.3GeV 2 in LO and NLO with PLUTO, L3, and ALEPH data as well as with the AFG results at Q 2 0 = 1GeV 2 in NLO systematic error. Finally, in this figure, our findings are compared with the LEP-DELPHI data [6] at scale Q 2 = 12 GeV 2 , in which the total error has been calculated from the quadratic sum of the statistical error and the systematic error. At the end of the presentation of the results for the photon structure function, in Fig. 8, we presented the results of the photon structure function at different initial scales Q 2 0 = 0.3 and 1 GeV 2 and compared them with the PLUTO [29], L3 [3], and ALEPH [7] data as well as with the AFG results [31] at the initial scale Q 2 0 = 1 GeV 2 . It is also noteworthy that the photon structure function increases by decreasing the initial scale, as shown in Fig. 8. For comparison, in Ref. [31], this increase at the peak of the figures is about 20%, but in our results, it is around 25%.
In Fig. 9, we plotted the photon distribution function obtained from Eq. (70) at LO and NLO in energy scales Q 2 = 5, 10, 20, 50, 100 and 1000 GeV 2 . Knowing that the photon distribution function has not been presented in the new data, we compared it with the s-quark and c-quark distribution functions of the photon at NLO provided by the GRV results. This distribution of the photon can be presented as a prediction of the DGLAP evolution equations.
In Tables 2 and 3, we presented our results, with notation "Our", of χ 2 -values, derived from our calculations and the experimental data and compared them with the resulting χ 2 -values of GRV and AFG. As can be seen in Table 2, our results at LO are in better agreement with experimental data than those obtained by GRV. However, at NLO, our results are only in agreement with the L3 and DELPHI data, which are again better than the GRV results achieved through the parametrization model. On the other hand, Table 3 shows the χ 2 -values of our results and the AFG results. Indeed, as shown, there is not much difference between these two results Fig. 9 The photon distribution function at scales Q 2 = 5, 10, 20, 50, 100, and 1000 GeV 2 in LO and NLO in comparison with the c-quark and s-quark distribution functions of the photon at NLO provided by the GRV results except for the experimental data released by the ALEPH collaboration. As illustrated, our analytical solutions for the photon structure function over a wide range of x and Q 2 values are in good agreement, especially at LO, with the experimental data provided by L3, OPAL, DELPHI, and PLUTO collaborations.
In conclusion, in this paper, we utilized the Laplace transform method for decoupling of the QCD⊗QED DGLAP evolution equations and extracted the singlet F  for the evolution. In the present paper, we also calculated the photon structure function F γ 2 (x, Q 2 ) using the Laplace transform method at the different values of Q 2 , derived from singlet F γ s (x, Q 2 ) and non-singlet F γ ns (x, Q 2 ) distribution functions. The solutions seem to be right, because the distribution functions of the photon and the photon structure function are in agreement with those from the available experimental data as well as with other parameterization models.

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: In this paper, we showed the results graphically and data are derived from the obtained equations.

The comparisons have just been shown numerically in several tables.]
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A
The Laplace integral is a powerful tool, formulated like the power series and Fourier series, to solve a wide variety of initial value problems. This integral was originally investigated in the pursuit of purely mathematical aims. The strategy is to transform the difficult differential equations into simple algebra problems (for a comprehensive discussion, see Ref. [35]). For a function F(x) defined on 0 ≤ x < ∞, the corresponding Laplace transform f (s) can be obtained by the following integral: where the parameter s may be real-valued or complex-valued and L is called the Laplace transform operator. The sufficient conditions for the existence of a Laplace transform of F(x) to f (s), for all s > a, are Furthermore, the inverse Laplace transform of a function f (s) is defined as follows: where L −1 is called the inverse Laplace transform operator and σ is a real number so that the contour path of the integration is in the region of convergence of f (s). The conditions for the existence of an inverse Laplace transform of f (s) to F(x) are as follows: (1) lim s→∞ f (s) = 0 (2) lim s→∞ s f (s) = finite.
The Mellin transform is a basic tool for analyzing the behavior of a function because of its scale invariance property. Indeed, this scale invariance property is analogous to the Fourier transform's shift-invariance property (for a comprehensive discussion, see Ref. [36]). The Mellin transform of an integrable function F(x) on (0, ∞) is defined as follows: where Also, the inverse Mellin transform is defined as where c ∈ (a, b) is any real number and M −1 is called the inverse Mellin transform operator. This integral exists when: (1) f (n) is regular in the infinite strip a < σ < b and, for any arbitrary small positive , f (n) → 0 in the strip a + ≤ σ ≤ b + as μ → ±∞. ∈ (a, b) is convergent and integrable.
In this paper, we used the following property in our calculations to simplify the equations: where H (

Appendix B
In Eqs.
In Eqs. where In Eqs. where Moreover, in Eqs.