A new numerical method for obtaining gluon distribution functions G(x, Q 2 )=xg(x, Q 2 ) , from the proton structure function F γp 2 (x, Q 2 )
aa r X i v : . [ h e p - ph ] J u l A new numerical method for obtaining gluon distribution functions G ( x, Q ) = xg ( x, Q ) , from the proton structure function F γp ( x, Q ) . Martin M. Block
Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208 (Dated: May 29, 2018)An exact expression for the leading-order (LO) gluon distribution function G ( x, Q ) = xg ( x, Q )from the DGLAP evolution equation for the proton structure function F γp ( x, Q ) for deep inelastic γ ∗ p scattering has recently been obtained [M. M. Block, L. Durand and D. W. McKay, Phys.Rev. D , 014031, (2009)] for massless quarks, using Laplace transformation techniques. Here, wedevelop a fast and accurate numerical inverse Laplace transformation algorithm, required to invertthe Laplace transforms needed to evaluate G ( x, Q ), and compare it to the exact solution. We obtainaccuracies of less than 1 part in 1000 over the entire x and Q spectrum. Since no analytic Laplaceinversion is possible for next-to-leading order (NLO) and higher orders, this numerical algorithmwill enable one to obtain accurate NLO (and NNLO) gluon distributions, using only experimentalmeasurements of F γp ( x, Q ). 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, accurateknowledge of gluon distribution functions at small Bjorken x will play a vital role in estimating backgrounds, andhence, 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 experi-mental data (mainly at small x ) on the proton structure function F γp ( x, Q ) measured in deep inelastic ep (or γ ∗ p )scattering, over a large domain of values of x and Q . The process starts with an initial Q , typically in the 1 to 2GeV range, and individual quark and gluon trial distributions parameterized as functions of x . The distributions areevolved to larger Q using the coupled integral-differential DGLAP equations [1, 2, 3], and the results used to predictthe measured quantities. The final distributions are then determined by adjusting the input parameters to obtain abest 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 ) does notappear in the experimentally accessible quantity F γp ( x, Q ), and is determined only through the quark distributionsin conjunction with the evolution equations. It is further not clear without detailed analysis [5, 8, 9, 10] how sensitivethe results are to the parameterizations of the initial parton distributions, or how well the gluon distribution is actuallydetermined.In a recent paper, Block, Durand and McKay (BDM) [11] used a Laplace transformation technique to obtaina leading order (LO) analytic gluon distribution function G ( x, Q ) = xg ( x, Q ) for massless quarks, directly froma global parameterization of the data on F γp ( x, Q ). The method uses only the LO DGLAP evolution equation[1, 2, 3] for F γp ( x, Q ) in the usual approximation in which the active quarks are treated as massless. In contrastto previous methods for determining G ( x, Q ), it does not require knowledge of the separate quark distributions inthe region in which structure function data exist, nor does it require the use of the evolution equation for G ( x, Q ),both considerable simplifications. In essence, the authors transform the LO DGLAP equation from Bjorken x -spaceto v -space, where v ≡ ln(1 /x ), and then Laplace transform the resulting equation. After solving for the Laplacetransform g ( s, Q ), 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 functionsneeded in the NLO DGLAP evolution equation [1, 2, 3] for F γp ( x, Q ) , the method can not be extended to NLOgluon 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 ( x, Q ) can be extendedto NLO (and NNLO) easily and accurately. Mathematically, Laplace inversion is an “ill-posed” problem and greatcare must be taken to insure its reliability for arbitrary Laplace transforms [12]. In this case, we can check ournumerical 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 ( x, Q ) for 4 massless quarks( u, d, s, c ) can be written as ∂F γp ( x, Q ) ∂ ln Q = α s π (cid:20) x Z x dzz F γp ( z, Q ) K qq (cid:16) xz (cid:17) + 209 x Z x dzz G ( z, Q ) K qg (cid:16) xz (cid:17)(cid:21) , (1)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 ( x, Q ) F γp ( x, Q ) ≡ ∂F γp ( x, Q ) ∂ ln Q − α s π x Z x dzz F γp ( z, Q ) K qq (cid:16) xz (cid:17) . (2)We finally write the DGLAP equation for the evolution of F γp ( x, Q ) as x Z x G ( z, Q ) K qg (cid:16) xz (cid:17) dzz = F ( x, Q ) , (3)where F ( x, Q ) ≡ (cid:16) α s π (cid:17) − F γp ( x, Q ) . (4)and the LO g → q splitting function is given by K qg ( x ) = 1 − x + 2 x . (5)At this point, BDM introduce the coordinate transformation v ≡ ln(1 /x ) , (6)and define functions ˆ G , ˆ K qg , and ˆ F in v -space byˆ G ( v, Q ) ≡ G ( e − v , Q )ˆ K qg ( v ) ≡ K qg ( e − v )ˆ F ( v, Q ) ≡ F ( e − v , Q ) . (7)Explicitly, from Eq. (5), we see that ˆ K qg ( v ) = 1 − e − v + 2 e − v . (8)Since ˆ F ( v, Q ) = Z v ˆ G ( w, Q ) e − ( v − w ) ˆ K qg ( v − w ) dw, (9)BDM note that in v -space, the DGLAP equation of Eq. (3) can be written as the convolution integralˆ F ( v, Q ) = Z v ˆ G ( w, Q ) ˆ H ( v − w ) dw, (10)where ˆ H ( v ) ≡ e − v ˆ K qg ( v )= e − v − e − v + 2 e − v . (11)The Laplace transform of ˆ H ( v ) is given by h ( s ), where h ( s ) ≡ L [ ˆ H ( v ); s ] = Z ∞ ˆ H ( v ) e − sv dv, with ˆ H ( v ) = 0 , v < . (12)The convolution theorem for Laplace transforms allows us to rewrite the right hand side of Eq. (10) as a product oftheir Laplace transforms g ( s ) and h ( s ), so that the Laplace transform of Eq. (10) is given by the algebraic equation f ( s, Q ) = g ( s, Q ) × h ( s ) . (13)Solving Eq. (13) for g , the Laplace transform of the gluon distribution function in s -space is given by g ( s, Q ) = f ( s, Q ) /h ( s ) . (14)In general, one is not able to calculate the inverse transform of g ( s, Q ) explicitly, if only because f ( s, Q ) is determinedby a numerical integral of the experimentally-determined function ˆ F ( v, Q ). However, regarding g ( s, Q ) as theproduct of the two functions f ( s, Q )and h − ( s ) and taking the inverse Laplace transform using the inverse of theconvolution theorem and the known inverse L − [ f ( s, Q ); v ] = ˆ F ( v, Q ), BDM find the analytic solutionˆ G ( v, Q ) = L − [ f ( s, Q ) × h − ( s ); v ] . (15)The calculation of h ( s ) and the inverse Laplace transform of h − ( s ) are straightforward, and the answer for LO,albeit singular, is given by BDM in terms of the Dirac delta function as L − [ h − ( s ); v ] = 3 δ ( v ) + δ ′ ( v ) − e − v/ √ " √ v + 2 cos " √ v . Thus, again using the convolution theorem, BDM find that the analytic solution for the LO gluon distribution ˆ G ( v, Q )for massless quarks isˆ G ( v, Q ) = 3 ˆ F ( v, Q ) + ∂ ˆ F ( v, Q ) ∂v − Z v ˆ F ( w, Q ) e − v − w ) / × √ " √
72 ( v − w ) + 2 cos " √
72 ( v − w ) dw. (16)The BDM solution depended critically upon the ability to find the analytic inverse Laplace transform of h − ( s ), whichwas 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 − is impossible—see Floratos et al. [13] for the NLO MS splitting function that is required. For this case, wemust be able to find a numerical inversion of the Laplace transform g ( s ) of the equivalent of Eq. (14) in order toobtain ˆ G ( v, Q ) and, ultimately, G ( x, Q ). The goal of this communication is to develop a suitable numerical Laplaceinversion algorithm. III. NUMERICAL INVERSION OF LAPLACE TRANSFORMS
In order to simplify our notation, we will now suppress the explicit dependence of ˆ G ( v, Q ) on Q , writing it asˆ G ( v ). If the Laplace transform of g ( s ) ≡ L [ ˆ G ( v ); s ] = R ∞ G ( v ) e − vs dv , where ˆ G ( v ) = 0 for v <
0, then the inverseLaplace transform is given by the complex Bromwich integralˆ G ( v ) ≡ L − [ g ( s ); v ] = 12 πi Z c + i ∞ c − i ∞ g ( s ) e vs ds, (17)where the real constant c is to the right of all singularities of g ( s ). We will further assume that we have made anappropriate coordinate translation in s so that c = 0, so that the equation is written asˆ G ( v ) ≡ L − [ g ( s ); v ] = 12 πi Z + i ∞− i ∞ g ( s ) e vs ds. (18)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 anew algorithm that takes advantage of very fast, arbitrarily high precision complex number arithmetic that is possibletoday in programs like Mathematica [14], making the inversion problem numerically tractable.First, we introduce a new complex variable z ≡ vs and rewrite Eq. (18) asˆ G ( v ) = 12 πiv Z + i ∞− i ∞ g (cid:16) zv (cid:17) e z dz. (19)We next make a rational approximation to e z , using the partial fraction expansion e z ≈ N X i =1 ω i z − α i , (20)which can be shown to have the following properties:1. Re α >
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, ω k z − α i + ¯ ω k z − ¯ α i , (21)is real for all real z .3. The expansion is identical to the Pad´e approximant with numerator equal to 2 N − N .4. The integrand vanishes faster than 1 /R as R → ∞ on the semi-circle of radius R that encloses the right handhalf 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 aclosed contour C by completing our integration path of the modified Bromwich integral in Eq. (19) with an infinitehalf circle to the right half of the complex plane. As mentioned earlier, g ( z/v ) has no singularities in this half of thecomplex plane. It is important to note that this contour is a clockwise path around the poles of Eq. (20), which comefrom our approximation to e z . What we need is the negative of it, i.e., the contour − C which is counterclockwise, sothat the poles are to our left as we traverse the contour − C . Therefore, we rewrite Eq. (19) asˆ G ( v ) ≈ πiv I C g (cid:16) zv (cid:17) N X i =1 ω i z − α i dz = − πiv N X i =1 I − C g (cid:16) zv (cid:17) ω i z − α i dz = − v N X i =1 Re [ ω i g ( α i /v )] . (22)To obtain Eq. (22), the final approximation to ˆ G ( v ), we used Cauchy’s theorem to equate the closed contour integralaround the path − C to 2 πi times the sum of the (complex) residues of the poles. Since the contour − C restrictsus to the right-hand half of the complex plane, no poles of g ( z/v ) were enclosed, but only the 2 N poles α i of theapproximation of e z . Using the properties cited above of the complex conjugate pairs—( ω i , α i ) and (¯ ω i , ¯ α i )—aftertaking only their real part and multiplying by 2, we have simultaneously insured that ˆ G ( v ) is real , yet only have hadto sum over half of the residues.Equation (22) has some very interesting properties:1. The 4 N coefficients ( α i , ω i ) are complex constants that are independent of v , only depending on the value of2 N used for the approximation, so that for a given 2 N , they only have to be evaluated once—in essence, theycan be tabulated and stored for later use.2. The Laplace transform of v n is given by n ! /s n +1 , where n is integer. Inserting G ( v ) = v n and g ( s ) = n ! /s n +1) into Eq. (22), we see that we have a set of 4 N equations, − n ! N X i (cid:18) ω i α n +1 i (cid:19) = 1 , n = 0 , , . . . , N − , (23)which also uniquely determine the 4 N complex constants ( α i , ω i ) , i = 1 , , . . . , N , although evaluating themdirectly from Eq. (23) is virtually impossible numerically, considering the ill-posed nature of these equations.The true power of Eq. (23) that it shows that the inversion algorithm of Eq. (22) is exact when G ( v ) is apolynomial of order ≤ N −
1, even though only N complex terms have to be evaluated in Eq. (22). This isreminiscent of the situation using Gauss-Legendre integration of order N , where there are 2 N constants, N Legendre zeroes and N weights, and the integration approximation is exact if the integrand is a polynomial oforder ≤ N − ω i alternate in sign and are exceedingly large—even for relatively modest 2 N ,making round-off a potentially serious problem. Thus, exceedingly high precision complex arithmetic is calledfor, often requiring 60 or more digits. However, this is not a serious problem—either in speed or complexity ofexecution—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 ) 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 havecalculated both the exact solution for ˆ G ( v, Q ) from Eq. (16) and the completely numerical approximation for ˆ G ( v )given in Eq. (22), using our inverse Laplace transformation algorithm given in Appendix A. For this purpose, we used N = 8 and p rec = 80 in the algorithm. The results for Q = 5 GeV and 100 GeV 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 ouralgorithm for the numerical inversion of Laplace transforms. As seen in both Figures, the agreement is striking overthe entire v range, which corresponds to the x interval 5 × − ≤ x ≤
1. At the highest v values (where numericalLaplace inversion approximations generally have the greatest problem), we find an accuracy of ∼ Q = 100 GeV . This error reflects both the numerical inaccuracy associated with the numerical integration term inthe exact solution ˆ G ( v, Q ) of Eq. (16), as well as the inaccuracies associated with the numerical Laplace inversionroutine 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 originalDGLAP equation from which we started, Eq. (3), i.e., x Z x G ( z, Q ) K qg (cid:16) xz (cid:17) dzz = F ( x, Q ) , (24)and numerically integrate its l.h.s., which depends on our numerical solution for ˆ G ( v, Q ) through G ( x, Q ) =ˆ G (ln(1 /x ) , Q ). We then compare it with the r.h.s., F ( x, Q ), which is independently known for all x and Q .This check becomes of primary importance when one does not have an analytical solution—the typical situation. Inthis case, it validated our conclusions about the numerical accuracy of our inversion procedure over the entire x and Q domain. In particular, for Q = 100 GeV , the ratio of the l.h.s. to the r.h.s. of Eq. (24) is unity to 1 part in5000 in the x range 3 × − < x < × − , and never rises to more than ∼ h − ( s ) in Eq. (16), theoverall behavior of ˆ G ( v, Q ), the Laplace inversion of f ( s, Q ) /h ( s ) in Eq. (16), is very smooth and therefore lendsitself readily to numerical inversion, as seen by the excellent agreement shown in Fig. 1 and Fig. 2 between the exactsolutions 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 becommented 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 v G L O H L H (cid:144) L
FIG. 1: LO gluon distribution functions ˆ G ( v, Q ) vs. v = ln(1 /x ), for virtuality Q = 5 GeV . The blue points are the exactLO solution of Eq. (16). The red curve is our numerical Laplace inversion solution. The agreement is excellent over the entire v range, which corresponds to the x -range, 5 × − ≤ x ≤ v G L O H L H (cid:144) L points are exact solution, curve is numerical solution
FIG. 2: LO gluon distribution functions ˆ G ( v, Q ) vs. v = ln(1 /x ), for virtuality Q = 100 GeV . The blue points are the exactLO solution of Eq. (16). The red curve is our numerical Laplace inversion solution. The agreement is excellent over the entire v range, which corresponds to the x -range, 5 × − ≤ x ≤ analytic solution using the same fit. Since NLO and NNLO DGLAP solutions are not amenable to having analyticsolutions—one can not analytically invert their h − ( s )—this technique will allow us to find very accurate numericalgluon solutions directly from experimental data, without having first to find quark distributions by solving coupledDGLAP equations. Further, we will be able to verify the numerical accuracy of our solutions by direct substitutioninto their generating equations. APPENDIX A:
MATHEMATICA
LAPLACE INVERSION ALGORITHM
NInverseLaplaceTransformBlock[g ,s ,v ,twoN ,prec ]:=Module[ { Omega,Alpha,M,p,den,r,num } ,(M=2*Ceiling[twoN/2];p=PadeApproximant[Exp[z], { z,0, { M-1,M }} ];den=Denominator[p];r=Roots[den==0,z]; Alpha=Table[r[[i,2]], { i,1,M } ];num=Numerator[p];hospital=num/D[den,z];Omega=SetPrecision[Table[hospital/.z->Alpha[[i]], { i,1,M,2 } ],prec+50];Alpha=SetPrecision[Table[Alpha[[i]], { i,1,M,2 } ],prec]);SetPrecision[-(2/v)Sum[Re[Omega[[i]] g/.s->Alpha[[i]]/v], { i,1,M/2 } ],30]] In the above algorithm, g = g ( s ), s = s , v = v , twoN =2 N in Eq. (22), and prec = precision of calculation. Typicalvalues are twoN = 8 and prec = 70. The algorithm, which is quite fast, returns the numerical value of ˆ G ( v ). It utilizesthat fact that the partial fraction approximation made for e z is equal to a Pad´e approximant whose numerator isorder 2 N − N .The algorithm first insures that M=twoN is even . It then constructs p , the Pad´e approximant whose numeratoris a polynomial of order M − M . It finds r , the complex roots of thedenominator, which are α i , the poles of Eq. (22). Using L’Hospital’s rule, it finds the residue ω i corresponding to thepole α i . At this point, all of the mathematics is symbolic. It next finds every other pair of ( α i , ω i ) to the desirednumerical accuracy; they come consecutively, i.e., α = ¯ α , ω = ¯ ω , α = ¯ α , ω = ¯ ω , etc. Finally, it takes thenecessary sums, again to the desired numerical accuracy, but only over half of the interval i = 1 , , . . . , N , by takingonly 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 achievearbitrarily high accuracy. If we define the accuracy as 1 − G ( v ) numerical /G ( v ) true , numerical tests on many differentfunctions shows that it goes to 0 for large as an inverse power law in . 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 —in the neighborhoodof 6 −
12, limiting the overall accuracy to be in the neighborhood of 10 − , which fortunately is ample for mostnumerical work. Typically, numerical integration routines are not accurate to better than ∼ − , and one can notuse ω ’s—which alternate in sign—that are larger than ∼ − , which occur for relatively small values of . Ofcourse, 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 g ( s ) = f ( s ) /h ( s ) . (B1)Although we can obtain an analytic solution for h ( s )—true for LO as well as NLO—we must find f ( s ) by numericalmeans. From Eq. (12), we see that the Laplace transform of ˆ F ( v ) is given by the integral f ( s ) = Z ∞ ˆ F ( v ) e − sv dv. (B2)Because of the ∞ at the upper limit of the integral, for numerical work it is useful to make the transformation v = − ln u , yielding h ( s ) = Z ˆ H ( − ln u ) u s − du. (B3)The upper limit of the integral is now finite, leading to the stable Mathematica algorithm:
NLaplaceTransform[H , vv , ss ] := NIntegrate[u^ (ss - 1)*H /.vv -> -Log[u], { u, 0, 1 } ,AccuracyGoal -> 15, MaxRecursion -> 15] APPENDIX C: GLOBAL PARAMETERIZATION OF F γp ( x, Q ) USING ZEUS STRUCTURE FUNCTIONDATA
Berger, Block and Tan [15] have parameterized the proton structure function F γp ( x, Q ) as F γp ( x, Q ) = (1 − x ) ( F P − x P + A ( Q ) ln (cid:20) x P x − x − x P (cid:21) + B ( Q ) ln (cid:20) x P x − x − x P (cid:21) ) . (C1)Here x P = 0 .
09 specifies the location in x of an approximate fixed point observed in the data where curves for different Q cross. At that point, ∂F γp ( x P , Q ) /∂ ln Q ≈ Q ; F P = F γp ( x P , Q ) = 0 .
41 is the common value of F γp .The Q dependence of F γp ( x, Q ) is given in those fits by A ( Q ) = a + a ln Q + a ln Q ,B ( Q ) = b + b ln Q + b ln Q . (C2)The fitted quantities and their errors are shown in Table I. TABLE I: Results of a 6-parameter fit to ZEUS F p ( x, Q ) structure function data [16] using the x and Q behaviors of Eq. (C1)and Eq. (C2), with Q in GeV . The renormalized χ per degree of freedom, taking into account the effects of the ∆ χ i max = 6cut [17], is given in the row labeled R × χ /d.f. The errors in the fitted parameters are multiplied by the appropriate r χ [17].Parameters Values a − . × − ± . × − a . × − ± . × − a . × − ± . × − b . × − ± . × − b . × − ± . × − b . × − ± . × − χ R × χ R × χ /d.f. 1.09 The fit to the data on F γp ( x, Q ) was restricted to the region x ≤ x P . In the absence of a global fit to the data inthis region, they simply extended their parameterization, piecewise, to the large- x region, using the form F γp ( x, Q ) = F P (cid:18) xx P (cid:19) µ ( Q ) (cid:18) − x − x P (cid:19) , x P < x ≤ , (C3)where the piecewise extension on the r.h.s. of Eq. (C3) is obviously continuous with the l.h.s. at x = x P , independentlyof µ ( Q ). The exponent µ ( Q ) is determined by requiring that the first derivatives with respect to x of the function inEqs. (C1) and the function on the r.h.s. of (C3) also match at x = x P . These results give the required parameterizationof F γp ( x, Q ) over all x , required in Eq. (4) to evaluate F ( x, Q ), and eventually, ˆ F ( v, Q ) of Eq. (7), needed to evaluatethe exact solution for ˆ G ( v, Q ) in Eq. (16), as well as calculating the Laplace transform f ( s, Q ) used in the numericalsolution of Eq. (14), which was calculated using the numerical algorithm of Appendix B. ACKNOWLEDGMENTS
Acknowledgments:
The author would like to thank Prof. L. Durand and Prof. Douglas W. McKay for theircontributions to portions of this work, and also to thank the Aspen Center for Physics for its hospitality during thetime parts of this work were done. [1] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. , 438 (1972).[2] G. Altarelli and G. Parisi, Nucl. Phys. B , 298 (1977).[3] Y. L. Dokshitzer, Sov. Phys. JETP , 641 (1977).[4] J. Pumplin et al. (CTEQ), J. High Energy Phys. , 012 (2002), hep-ph/0201195.[5] W. K. Tung, H. L. Lai, A. Belyaev, J. Pumplin, D. Stump, and C.-P. Yuan, J. High Energy Phys. , 053 (2007),hep-ph/0611254.[6] A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S. Thorne, Eur. Phys. J. C , 73 (2002), hep-ph/0110215.[7] A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S. Thorne, Phys. Lett. B , 61 (2004), hep-ph/0410230.[8] A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S. Thorne, Eur. Phys. J. C , 325 (2004), hep-ph/0308087.[9] D. Stump et al. (ZEUS), Phys. Rev. D , 014012 (2002), hep-ph/0101051.[10] D. Stump et al. (ZEUS), Phys. Rev. D , 014013 (2002), hep-ph/0101032.[11] M. M. Block, L. Durand, and D. W. McKay, Phys. Rev. D , 014031 (2008), arXiv:0808.0201 [hep-ph].[12] U. Graf, Applied Laplace Transforms and z-Transforms for Scidentists and Engineers, Birkhause Verlag, Basel, Switzerland.“This is not a so-called well-posed problem. The inverse Laplace transform is an unbounded operator which gives rise toa numerically ill-conditioned problem”, p. ix (2003).[13] E. G. Floratos, C. Kounas, and R. Lacaze, Nucl. Phys. B , 417 (1981).[14] Mathematica , 242001 (2007), hep-ph/0703003.[16] S. Chekanov et al. (ZEUS), Eur. Phys. J. C , 443 (2001).[17] M. M. Block, Nucl. Inst. and Meth. A.556