A new numerical method for obtaining gluon distribution functions $G(x,Q^2)=xg(x,Q^2)$, from the proton structure function $F_2^{\gamma p}(x,Q^2)$

An exact expression for the leading-order (LO) gluon distribution function $G(x,Q^2)=xg(x,Q^2)$ from the DGLAP evolution equation for the proton structure function $F_2^{\gamma p}(x,Q^2)$ for deep inelastic $\gamma^* p$ scattering has recently been obtained [M. M. Block, L. Durand and D. W. McKay, Phys. Rev. D{\bf 79}, 014031, (2009)] for massless quarks, using Laplace transformation techniques. Here, we develop a fast and accurate numerical inverse Laplace transformation algorithm, required to invert the Laplace transforms needed to evaluate $G(x,Q^2)$, and compare it to the exact solution. We obtain accuracies of less than 1 part in 1000 over the entire $x$ and $Q^2$ spectrum. Since no analytic Laplace inversion is possible for next-to-leading order (NLO) and higher orders, this numerical algorithm will enable one to obtain accurate NLO (and NNLO) gluon distributions, using only experimental measurements of $F_2^{\gamma p}(x,Q^2)$.


I. INTRODUCTION
The quark and gluon distributions in hadrons play a key role in our understanding of Standard Model processes, in our predictions for such processes at accelerators, and in our searches for new physics. In particular, accurate knowledge of gluon distribution functions at small Bjorken x will play a vital role in estimating backgrounds, and hence, our ability to search for new physics at the Large Hadron Collider.
The gluon and quark distribution functions have traditionally been determined simultaneously by fitting experimental data (mainly at small x) on the proton structure function F γp 2 (x, Q 2 ) measured in deep inelastic ep (or γ * p) scattering, over a large domain of values of x and Q 2 . The process starts with an initial Q 2 0 , typically in the 1 to 2 GeV 2 range, and individual quark and gluon trial distributions parameterized as functions of x. The distributions are evolved to larger Q 2 using the coupled integral-differential DGLAP equations [1,2,3], and the results used to predict the measured quantities. The final distributions are then determined by adjusting the input parameters to obtain a best fit to the data. For recent determinations of the gluon and quark distributions, see [4,5,6,7].
This procedure is rather indirect, especially so in the case of the gluon: the gluon distribution G(x, Q 2 ) does not appear in the experimentally accessible quantity F γp 2 (x, Q 2 ), and is determined only through the quark distributions in conjunction with the evolution equations. It is further not clear without detailed analysis [5,8,9,10] how sensitive the results are to the parameterizations of the initial parton distributions, or how well the gluon distribution is actually determined.
In a recent paper, Block, Durand and McKay (BDM) [11] used a Laplace transformation technique to obtain a leading order (LO) analytic gluon distribution function G(x, Q 2 ) = xg(x, Q 2 ) for massless quarks, directly from a global parameterization of the data on F γp 2 (x, Q 2 ). The method uses only the LO DGLAP evolution equation [1,2,3] for F γp 2 (x, Q 2 ) in the usual approximation in which the active quarks are treated as massless. In contrast to previous methods for determining G(x, Q 2 ), it does not require knowledge of the separate quark distributions in the region in which structure function data exist, nor does it require the use of the evolution equation for G(x, Q 2 ), both considerable simplifications. In essence, the authors transform the LO DGLAP equation from Bjorken x-space to v-space, where v ≡ ln(1/x), and then Laplace transform the resulting equation. After solving for the Laplace transform g(s, Q 2 ), they are able to analytically invert the Laplace transform back into v-space, and eventually, to x-space. Unfortunately, because of the considerable complexities of next-to-leading order (NLO) splitting functions needed in the NLO DGLAP evolution equation [1,2,3] for F γp 2 (x, Q 2 ) , the method can not be extended to NLO gluon distributions, because of the impossibility of analytically inverting the required Laplace transform.
The purpose of this note is to derive a new and fast algorithm for accurate numerical inversion of Laplace transforms, so that the BDM method [11] for direct evaluation of gluon distributions from knowledge of F γp 2 (x, Q 2 ) can be extended to NLO (and NNLO) easily and accurately. Mathematically, Laplace inversion is an "ill-posed" problem and great care must be taken to insure its reliability for arbitrary Laplace transforms [12]. In this case, we can check our numerical Laplace inversion routine directly by comparing its results to the exact LO gluon solution of Ref. [11].

II. THE EXACT LO SOLUTION
The LO DGLAP equation for the evolution of the proton structure function F γp 2 (x, Q 2 ) for 4 massless quarks (u, d, s, c) can be written as where K qq (x) and K qg (x) are the LO splitting functions and α s is the strong running coupling constant. Following BDM [11], we introduce F γp We finally write the DGLAP equation for the evolution of F γp 2 (x, Q 2 ) as where and the LO g → q splitting function is given by At this point, BDM introduce the coordinate transformation and define functionsĜ,K qg , andF in v-space bŷ Explicitly, from Eq. (5), we see thatK SinceF BDM note that in v-space, the DGLAP equation of Eq. (3) can be written as the convolution integral whereĤ The Laplace transform ofĤ(v) is given by h(s), where The convolution theorem for Laplace transforms allows us to rewrite the right hand side of Eq. (10) as a product of their Laplace transforms g(s) and h(s), so that the Laplace transform of Eq. (10) is given by the algebraic equation Solving Eq. (13) for g, the Laplace transform of the gluon distribution function in s-space is given by In general, one is not able to calculate the inverse transform of g(s, Q 2 ) explicitly, if only because f (s, Q 2 ) is determined by a numerical integral of the experimentally-determined functionF (v, Q 2 ). However, regarding g(s, Q 2 ) as the product of the two functions f (s, Q 2 )and h −1 (s) and taking the inverse Laplace transform using the inverse of the convolution theorem and the known inverse The calculation of h(s) and the inverse Laplace transform of h −1 (s) are straightforward, and the answer for LO, albeit singular, is given by BDM in terms of the Dirac delta function as Thus, again using the convolution theorem, BDM find that the analytic solution for the LO gluon distributionĜ(v, Q 2 ) for massless quarks iŝ The BDM solution depended critically upon the ability to find the analytic inverse Laplace transform of h −1 (s), which was possible for the LO case of the splitting function K qg (x). Unfortunately, for NLO or higher, the splitting function is so complicated that an analytic inversion for the NLO h −1 is impossible-see Floratos et al. [13] for the NLO MS splitting function that is required. For this case, we must be able to find a numerical inversion of the Laplace transform g(s) of the equivalent of Eq. (14) in order to obtainĜ(v, Q 2 ) and, ultimately, G(x, Q 2 ). The goal of this communication is to develop a suitable numerical Laplace inversion algorithm.

III. NUMERICAL INVERSION OF LAPLACE TRANSFORMS
In order to simplify our notation, we will now suppress the explicit dependence ofĜ(v,

then the inverse Laplace transform is given by the complex Bromwich integral
where the real constant c is to the right of all singularities of g(s). We will further assume that we have made an appropriate coordinate translation in s so that c = 0, so that the equation is written aŝ Our goal is to numerically solve Eq. (18). The inverse Laplace transform is essentially determined by the behavior of g(s) near its singularities, and thus is an ill-conditioned or ill-posed numerical problem. We suggest in this note a new algorithm that takes advantage of very fast, arbitrarily high precision complex number arithmetic that is possible today in programs like Mathematica [14], making the inversion problem numerically tractable.
First, we introduce a new complex variable z ≡ vs and rewrite Eq. (18) aŝ We next make a rational approximation to e z , using the partial fraction expansion which can be shown to have the following properties: 1. Re α 1 > 0, so that its poles are all in the right-hand half of the complex plane.
2. N distinct complex conjugate pairs of the complex numbers (ω i , α i ), such that the sum of the k th pair, is real for all real z.
3. The expansion is identical to the Padé approximant with numerator equal to 2N − 1 and denominator equal to 2N .
4. The integrand vanishes faster than 1/R as R → ∞ on the semi-circle of radius R that encloses the right hand half of the complex plane, since the approximation vanishes as 1/R and g(s) that corresponds to a non-singular G(v) also vanishes for R → ∞.
Since g(z/v) must vanish for z → ∞ and our approximation for e z in Eq. (20) vanishes for z → ∞, we can form a closed contour C by completing our integration path of the modified Bromwich integral in Eq. (19) with an infinite half circle to the right half of the complex plane. As mentioned earlier, g(z/v) has no singularities in this half of the complex plane. It is important to note that this contour is a clockwise path around the poles of Eq. (20), which come from our approximation to e z . What we need is the negative of it, i.e., the contour −C which is counterclockwise, so that the poles are to our left as we traverse the contour −C. Therefore, we rewrite Eq. (19) aŝ To obtain Eq. (22), the final approximation toĜ(v), we used Cauchy's theorem to equate the closed contour integral around the path −C to 2πi times the sum of the (complex) residues of the poles. Since the contour −C restricts us to the right-hand half of the complex plane, no poles of g(z/v) were enclosed, but only the 2N poles α i of the approximation of e z . Using the properties cited above of the complex conjugate pairs-(ω i , α i ) and (ω i ,ᾱ i )-after taking only their real part and multiplying by 2, we have simultaneously insured thatĜ(v) is real , yet only have had to sum over half of the residues. Equation (22) has some very interesting properties: 1. The 4N coefficients (α i , ω i ) are complex constants that are independent of v, only depending on the value of 2N used for the approximation, so that for a given 2N , they only have to be evaluated once-in essence, they can be tabulated and stored for later use.
The true power of Eq. (23) that it shows that the inversion algorithm of Eq. (22) is exact when G(v) is a polynomial of order ≤ 4N − 1, even though only N complex terms have to be evaluated in Eq. (22). This is reminiscent of the situation using Gauss-Legendre integration of order N , where there are 2N constants, N Legendre zeroes and N weights, and the integration approximation is exact if the integrand is a polynomial of order ≤ 2N − 1.
3. The real parts of the residues ω i alternate in sign and are exceedingly large-even for relatively modest 2N , making round-off a potentially serious problem. Thus, exceedingly high precision complex arithmetic is called for, often requiring 60 or more digits. However, this is not a serious problem-either in speed or complexity of execution-for an algorithm written in a program such as Mathematica [14].
A concise inversion algorithm in Mathematica that implements Eq. (22) is given in Appendix A and a one line Mathematica algorithm for finding a numerical Laplace transform is given in Appendix B.
We will now test the accuracy of our numerical Laplace inversion algorithm by comparing its results with Eq. (16), the exact LO G(v, Q 2 ) of BDM.

IV. COMPARISON OF EXACT SOLUTION AND NUMERICAL LAPLACE INVERSION RESULTS
Using the Berger, Block and Tan [15] fit to the experimental ZEUS [16] data shown in Appendix C, we have calculated both the exact solution forĜ(v, Q 2 ) from Eq. (16) and the completely numerical approximation forĜ(v) given in Eq. (22), using our inverse Laplace transformation algorithm given in Appendix A. For this purpose, we used 2N = 8 and prec = 80 in the algorithm. The results for Q 2 = 5 GeV 2 and 100 GeV 2 are shown in Fig. 1 and Fig.  2, respectively. The blue points are the exact solution and the red curves are the numerical solution that uses our algorithm for the numerical inversion of Laplace transforms. As seen in both Figures, the agreement is striking over the entire v range, which corresponds to the x interval 5 × 10 −7 ≤ x ≤ 1. At the highest v values (where numerical Laplace inversion approximations generally have the greatest problem), we find an accuracy of ∼ 1 part in 5000, for Q 2 = 100 GeV 2 . This error reflects both the numerical inaccuracy associated with the numerical integration term in the exact solutionĜ(v, Q 2 ) of Eq. (16), as well as the inaccuracies associated with the numerical Laplace inversion routine of Appendix A and the numerical Laplace transformation routine of Appendix B.
Another independent method of checking the numerical accuracv of the entire procedure is to go back to the original DGLAP equation from which we started, Eq. (3), i.e., and numerically integrate its l.h.s., which depends on our numerical solution forĜ(v, Q 2 ) through G(x, Q 2 ) = G(ln(1/x), Q 2 ). We then compare it with the r.h.s., F (x, Q 2 ), which is independently known for all x and Q 2 . This check becomes of primary importance when one does not have an analytical solution-the typical situation. In this case, it validated our conclusions about the numerical accuracy of our inversion procedure over the entire x and Q 2 domain. In particular, for Q 2 = 100 GeV 2 , the ratio of the l.h.s. to the r.h.s. of Eq. (24) is unity to 1 part in 5000 in the x range 3 × 10 −4 < x < 3 × 10 −2 , and never rises to more than ∼ 1 part in 200 outside this range; these excursions from unity include the errors generated in the numerical integration of the l.h.s of Eq. (24). In contrast to the highly singular behavior of the analytic Laplace inverse of inversion of h −1 (s) in Eq. (16), the overall behavior ofĜ(v, Q 2 ), the Laplace inversion of f (s, Q 2 )/h(s) in Eq. (16), is very smooth and therefore lends itself readily to numerical inversion, as seen by the excellent agreement shown in Fig. 1 and Fig. 2 between the exact solutions and the numerical inverse Laplace transforms.
Although our inversion routine was specifically developed to invert Laplace transforms of gluon distributions, it clearly has a wide variety of applications in solving both integral and differential equations, which will not be commented on further in this communication.

V. CONCLUSIONS
We have achieved high numerical accuracy in obtaining LO gluon distributions from fits to experimental data, using a numerical Laplace inversion routine developed for this purpose, and have compared it successfully to the exact analytic solution using the same fit. Since NLO and NNLO DGLAP solutions are not amenable to having analytic solutions-one can not analytically invert their h −1 (s)-this technique will allow us to find very accurate numerical gluon solutions directly from experimental data, without having first to find quark distributions by solving coupled DGLAP equations. Further, we will be able to verify the numerical accuracy of our solutions by direct substitution into their generating equations. In the above algorithm, g = g(s), s = s, v = v, twoN=2N in Eq. (22), and prec = precision of calculation. Typical values are twoN = 8 and prec = 70. The algorithm, which is quite fast, returns the numerical value ofĜ(v). It utilizes that fact that the partial fraction approximation made for e z is equal to a Padé approximant whose numerator is order 2N − 1 and whose denominator is order 2N .
The algorithm first insures that M=twoN is even. It then constructs p, the Padé approximant whose numerator is a polynomial of order M − 1 and denominator a polynomial of order M . It finds r, the complex roots of the denominator, which are α i , the poles of Eq. (22). Using L'Hospital's rule, it finds the residue ω i corresponding to the pole α i . At this point, all of the mathematics is symbolic. It next finds every other pair of (α i , ω i ) to the desired numerical accuracy; they come consecutively, i.e., α 1 =ᾱ 2 , ω 1 =ω 2 , α 3 =ᾱ 4 , ω 3 =ω 4 , etc. Finally, it takes the necessary sums, again to the desired numerical accuracy, but only over half of the interval i = 1, 3, . . . , N , by taking only the real part and multiplying by 2.
If g(s), the input to the algorithm, is an analytic relation and v is a pure number (from the point of view of Mathematica, 31/10 is a pure number, but 3.1 is not), then, for sufficiently high values of prec, you can achieve arbitrarily high accuracy. If we define the accuracy as 1 − G(v) numerical /G(v) true , numerical tests on many different functions shows that it goes to 0 for large 2N as an inverse power law in 2N. On the other hand, if g(s) is obtained numerically, either from having to find the Laplace transform f (s) and/or h(s) from numerical integration techniques, the overall accuracy of inversion is limited by the need to only use relatively small values of 2N-in the neighborhood of 6 − 12, limiting the overall accuracy to be in the neighborhood of 10 −5 , which fortunately is ample for most numerical work. Typically, numerical integration routines are not accurate to better than ∼ 10 −5 , and one can not use ω's-which alternate in sign-that are larger than ∼ 10 14 − 10 16 , which occur for relatively small values of 2N. Of course, this is not a limitation if g(s) is able to be expressed in closed form.

APPENDIX B: MATHEMATICA NUMERICAL LAPLACE TRANSFORM ALGORITHM
We note from Eq. (14) that the function that we must invert is Although we can obtain an analytic solution for h(s)-true for LO as well as NLO-we must find f (s) by numerical means. From Eq. (12), we see that the Laplace transform ofF (v) is given by the integral Because of the ∞ at the upper limit of the integral, for numerical work it is useful to make the transformation v = − ln u, yielding APPENDIX C: GLOBAL PARAMETERIZATION OF F γp 2 (x, Q 2 ) USING ZEUS STRUCTURE FUNCTION DATA Berger, Block and Tan [15] have parameterized the proton structure function F γp 2 (x, Q 2 ) as Here x P = 0.09 specifies the location in x of an approximate fixed point observed in the data where curves for different Q 2 cross. At that point, ∂F γp 2 (x P , Q 2 )/∂ ln Q 2 ≈ 0 for all Q 2 ; F P = F γp 2 (x P , Q 2 ) = 0.41 is the common value of F γp 2 . The Q 2 dependence of F γp 2 (x, Q 2 ) is given in those fits by A(Q 2 ) = a 0 + a 1 ln Q 2 + a 2 ln 2 Q 2 , The fitted quantities and their errors are shown in Table I.  [16] using the x and Q 2 behaviors of Eq. (C1) and Eq. (C2), with Q 2 in GeV 2 . The renormalized χ 2 min per degree of freedom, taking into account the effects of the ∆χ 2 i max = 6 cut [17], is given in the row labeled R × χ 2 min /d.f. The errors in the fitted parameters are multiplied by the appropriate rχ2 [17]. The fit to the data on F γp 2 (x, Q 2 ) was restricted to the region x ≤ x P . In the absence of a global fit to the data in this region, they simply extended their parameterization, piecewise, to the large-x region, using the form where the piecewise extension on the r.h.s. of Eq. (C3) is obviously continuous with the l.h.s. at x = x P , independently of µ(Q 2 ). The exponent µ(Q 2 ) is determined by requiring that the first derivatives with respect to x of the function in Eqs. (C1) and the function on the r.h.s. of (C3) also match at x = x P . These results give the required parameterization of F γp 2 (x, Q 2 ) over all x, required in Eq. (4) to evaluate F (x, Q 2 ), and eventually,F (v, Q 2 ) of Eq. (7), needed to evaluate the exact solution forĜ(v, Q 2 ) in Eq. (16), as well as calculating the Laplace transform f (s, Q 2 ) used in the numerical solution of Eq. (14), which was calculated using the numerical algorithm of Appendix B.