{'fraction_non_alphanumeric': 0.05296318172868675, 'fraction_numerical': 0.018655247095885015, 'mean_word_length': 4.042690493919086, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decoupling the NLO coupled DGLAP evolution equations: an analytic solution to pQCD

Using repeated Laplace transform techniques, along with newly-developed accurate numerical inverse Laplace transform algorithms, we transform the coupled, integral-differential NLO singlet DGLAP equations first into coupled differential equations, then into coupled algebraic equations, which we can solve iteratively. After Laplace inverting the algebraic solution analytically, we numerically invert the solutions of the decoupled differential equations. Finally, we arrive at the decoupled NLO evolved solutions F_s(x,Q^2)=calF_s(F_{s0}(x),G_0(x)) and G(x,Q^2)=calG(F_{s0}(x),G_0(x)), where calF_s and calG are known functions - determined using the DGLAP splitting functions up to NLO in the strong coupling constant alpha_s(Q^2). The functions F_{s0}(x)=F_s(x,Q_0^2) and G_0(x)=G(x,Q_0^2) are the starting functions for the evolution at Q_0^2. This approach furnishes us with a new tool for readily obtaining, independently, the effects of the starting functions on either the evolved gluon or singlet structure functions, as a function of both Q^2 and Q_0^2. It is not necessary to evolve coupled integral-differential equations numerically on a two-dimensional grid, as is currently done. The same approach can be used for NLO non-singlet distributions where it is simpler, only requiring one Laplace transform. We make successful NLO numerical comparisons to two non-singlet distributions, using NLO quark distributions published by the MSTW collaboration, over a large range of x and Q^2. Our method is readily generalized to higher orders in the strong coupling constant alpha_s(Q^2).

where Fs and G are known functions-determined using the DGLAP splitting functions up to NLO in the strong coupling constant αs(Q 2 ). The functions Fs0(x) ≡ Fs(x, Q 2 0 ) and G0(x) ≡ G(x, Q 2 0 ) are the starting functions for the evolution at Q 2 = Q 2 0 . This approach furnishes us with a new tool for readily obtaining, independently, the effects of the starting functions on either the evolved gluon or singlet structure functions, as a function of both Q 2 and Q 2 0 . It is not necessary to evolve coupled integral-differential equations numerically on a two-dimensional grid, as is currently done. The same approach can be used for NLO non-singlet distributions where it is simpler, only requiring one Laplace transform. We make successful NLO numerical comparisons to two non-singlet distributions, using NLO quark distributions published by the MSTW collaboration [3], over a large range of x and Q 2 . Our method is readily generalized to higher orders in the strong coupling constant αs(Q 2 ).

I. INTRODUCTION
In order to interpret the experimental results at the Large Hadron Collider-in the search for new physicsaccurate knowledge of gluon distribution functions at small Bjorken x and large virtuality Q 2 plays a vital role in estimating QCD backgrounds and in calculating gluon-initiated processes. The gluon and quark distribution functions have traditionally been determined simultaneously by fitting experimental data on neutral-and charged-current deep inelastic scattering processes and some jet data over a large domain of values of x and Q 2 . The distributions at small x and large Q 2 are determined mainly by the proton structure function F γp 2 (x, Q 2 ) measured in deep inelastic ep (or γ * p) scattering. The fitting process starts with an initial Q 2 0 , typically less than m 2 c , the square of the c quark mass of ≈ 2 GeV 2 , and individual quark and gluon initial distributions which are parameterized with pre-determined shapes in x determined by a set of adjustable input parameters-given as functions of x for the chosen Q 2 0 . The distributions are then evolved numerically on a two-dimensional grid in x and Q 2 to larger Q 2 using the coupled integral-differential DGLAP equations [4][5][6], typically in leading order (LO) and next-to-leading order (NLO), 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. This procedure is very indirect in the case of the gluon: the gluon distribution G(x, Q 2 ) = xg(x, Q 2 ) does not contribute directly to the accurately determined structure function F γp 2 (x, Q 2 ), and is determined only through the quark distributions in conjunction with the evolution equations, or at large x, from jet data. For recent determinations of the gluon and quark distributions, see [3,[7][8][9][10].
In the following, we will summarize our analytic method that determines the singlet structure function F s (x, Q 2 ) and G(x, Q 2 ) directly and individually, using as input F s0 (x) ≡ F s (x, Q 2 0 ) and G 0 (x) ≡ G(x, Q 2 0 ), where Q 2 0 is arbitrary, with the guarantee that each distribution individually satisfies the NLO coupled DGLAP equations.
The method is extended to calculate NLO non-singlet functions, so that we can also find individual quark and gluon distributions analytically in terms of the starting distributions of the individual quark and gluon distributions.
We will also give some numerical examples for non-singlet NLO valence quark distributions, comparing them to the MSTW [3] published NLO valence quark distributions.

II. NLO SINGLET SECTOR
Our approach uses an unusual application of multiple Laplace transforms [11,12]. In this note, we use double Laplace transforms, first transforming the coupled DGLAP integral-differential equations into a set of coupled differential equations in Laplace space, and finally, into a set of coupled algebraic equations in a second Laplace space. We then solve the coupled algebraic equations in this second Laplace space. To obtain our final results, we must invert the Laplace transforms. The second transform to the algebraic equation space is analytically invertible to the space in which we had the coupled differential equations. The final inversion, from this Laplace space back to our initial space, must be obtained by numerical inverse Laplace transformations [1,2].
We first introduce the variable v ≡ ln(1/x) into the NLO coupled DGLAP equations. This turns them into coupled convolution equations in v space, which, after introducing a new variable τ (Q 2 , , are readily Laplace transformed to obtain a set of coupled homogeneous first-order differential equations in the variable τ . The parameters of these transformed equations are known functions of s, the Laplace-space variable. These equations are then Laplace transformed a second time, essentially transforming the variable τ of the coupled differential equations into a new Laplace variable U , with the resulting equations being coupled algebraic equations in U -with s again being a parameter-which are then solved iteratively. These solutions, in s and U , are analytically Laplace inverted back to variables s and τ . Using fast and accurate numerical inverse Laplace transform algorithms [1,2], we transform the solutions back into v space, and, finally, into Bjorken x-space to obtain F s (x, Q 2 ) = F s (F s0 (x), G 0 (x)) and G(x, Q 2 ) = G(F s0 (x), G 0 (x)), where the functions F and G are determined by the splitting functions in the DGLAP equations.
A similar method was used in an earlier paper [13] in which we obtained the decoupled solutions in LO for both the singlet and the non-singlet sector, using only one Laplace transform. The τ dependence in that case was trivial, and the decoupled equations could be solved directly. The extra Laplace transform that appears in the present work is necessitated by the nontrivial dependence of the NLO terms on τ .
Our method is readily generalized to all orders in the strong coupling constant, but for brevity we limit ourselves to NLO in this paper. We write the coupled NLO DGLAP equations [11,12] schematically, using the convolution symbol ⊗, as The P 0 qq (x), P 0 qg (x), P 0 gq (x) and P 0 gg (x) used in Eq. (1) and Eq. (2) are the LO singlet splitting function and the P 1 qq (x), P 1 qg (x), P 1 gq (x) and P 1 gg (x) are the NLO singlet splitting functions, with α s (Q 2 ) the NLO running strong coupling constant. It is standard procedure to construct α s (Q 2 ) assuming three massless quarks, u, d and s, below the c-quark threshold, adjusting the QCD parameter Λ at each successive threshold in Q 2 , i,e, at Q 2 = M 2 c and M 2 b , so that α s (Q 2 ) remains continuous when the number of quarks changes as heavy c and b quarks begin to contribute.
Introducing the variable changes and the notationF we rewrite the above DGLAP equations in terms of the convolution integrals whereF The splitting functions P qq (e −v ) and P gg (e −v ) in the integrals Eq. (5) and Eq. (6) involve the distribution 1/(1 − e −v ) + . The integrals involving this term can be transformed to expressions that involve derivatives ofF s orĜ without the appearance of the singular factor 1/(1 − e −v ), for example, to integrals of the form After this change, all of the integrals in Eqs. (5) and (6) are normal convolutions. By making a Laplace transform in v, we can factor these integrals, since the Laplace transform of a convolution is the product of the Laplace transform of the factors, i.e., Defining the Laplace transforms and noting that we can factor the Laplace transforms of Eq. (5) and Eq. (6) into two coupled ordinary first order differential equations in the variable τ in Laplace space s with τ -dependent coefficients. These can be written as where we recall that the the Q 2 dependence is through the function τ , i.e., τ (Q 2 , The LO coefficients Φ LO and Θ LO are given by [13] Φ Here ψ(x) in Eqs. (15) and (17) is the digamma function and γ E = 0.5772156 . . . is Euler's constant, quantities that are introduced in the Laplace transform of the LO terms involving the distribution 1/(1 − e −v ) + discussed above. The evaluation of the NLO coefficients is straightforward, but too lengthy to be included in this note, and will be given, in the future, when we make numerical evaluations of F s (x, Q 2 ) and G(x, Q 2 ) in NLO.
In the case of LO, the τ dependence of the equations is trivial, and the equations can be solved simply [13], as already noted. The extra explicit dependence of the NLO terms on the right-hands of these equations on τ prevents a similar construction here. In order to decouple and solve Eq. (13) and Eq. (14), we Laplace transform them a second time-this time with respect to the variable τ -into U space, i.e., we let where now s is simply a parameter in U space.
In U space, we now write the final desired coupled algebraic equations for F (s, U ) and G(s, U ) as For brevity, we replace the NLO αs(τ ) 4π by a(τ ) in Eq. (21) and Eq. (22). We can numerically show that an excellent approximation to a(τ ) ≡ αs(τ ) 4π , accurate to a few parts in 10 4 , is given by the expression where the constants a 0 , a 1 , b 1 are found by a least squares fit to a(τ ). We note in passing that this approximation is inspired by the fact that in LO, α s,LO (τ ) is exactly given by the form α s,LO (Q 2 0 )e −bτ . Using the value of a(τ ) given by Eq.
where b 0 = 0 is understood in Eq. (23) and Eq. (24). After introducing the simplifying notation we can finally rewrite Eq. (21) and Eq. (22) as which we will solve iteratively, using a 1 as an expansion parameter. We note that the Φ's and Θ's, as defined above, contain both LO and NLO terms. We further point out that Eq. (27) and Eq. (28) are completely symmetric under the simultaneous transformations f ↔ g and F ↔ G. We finally remark that a 1 , the NLO expansion parameter in our iterative solution of Eqs. We next consider the simple solutions to Eq. (27) and Eq. (28), called F 1 (s, U ) and G 1 (s, U ), that result from setting a 1 = 0, i.e., the equations whose solutions are The denominator D(U, s) in Eqs. (31) and (32) is just the determinant of the coefficients of F (s, U ) and G(s, U ) in Eqs. (27) and (28), where . The zeros of D(U, s) lead to simple poles in F 1 and G 1 in the U plane. These functions have no other singularities, and decrease as 1/|U | for |U | → ∞. The inverse Laplace transforms of F 1 (s, U ) and G 1 (s, U ), denoted by f 1 (s, τ ) and g 1 (s, τ ), are therefore well defined and simple to calculate. We will write them as where the coefficient functions in the solution are We comment that this solution has small NLO terms in it, arising from the a 0 term in Eq. (25) and Eq. (26). However, it is identical in form to the LO solution we gave in Ref. [13], and reduces to it if we set a 0 = 0. We next construct an iterative solution to Eqs. (27) and (28) for F and G. We start by substituting the known functions F 1 and G 1 for F and G on the right hand sides of the equations, and then re-solve the equations to obtain the next approximations F 2 and G 2 for F and G, and then repeat the process. For the first step, we make the replacements on the right-hand sides of Eqs. (27) and (28) to obtain our first iterative equations for F 2 (s, U ) and G 2 (s, U ), The functions F 1 (s, U ) and G 1 (s, U ) are given analytically by Eq. (31) and Eq. (32), respectively, so that we know them at the argument (s, U + b 1 ), needed in the right hand sides of our iterative equations.
Since the functions on the right-hand sides of Eqs. (40) and (41) are known, and the left-hand sides have the same structure as Eqs. (29) and (30), their solutions can be obtained by the substitutions on the right hand sides of Eqs. (31) and (32). The leading f 0 (s) and g 0 (s) reproduce F 1 (s, U ) and G 1 (s, U ). The added terms, proportional to the expansion parameter a 1 , are more complicated expressions that are rational functions in U , whose numerators are a second-order polynomial and whose denominators are the factorable product D(U, s)D(U + b 1 , s). The functions F 2 and G 2 therefore have an extra pole in U that is displaced along the real axis from the poles of F 1 and G 1 by the amount b 1 . Since this is the only new singularity, whose terms decrease at least as rapidly as 1/U 2 for U → ∞, the overall behavior of the iterated solution decreases at least as rapidly as 1/|U |, so that the inverse Laplace transforms needed can be calculated analytically. We can again write the results for the inverse transforms f 2 (s, τ ) and g 2 (s, τ ) in terms of the initial distributions f 0 (s) and g 0 (s) as in Eq. (34), but with the coefficient functions k now sums of the original expressions in Eqs. (35)-(38) and terms that depend linearly on the coefficient a 1 in Eq. (23) as well as on s and τ . The coefficient a 0 has been incorporated into the definitions of the Φ's and Θ's in Eqs. (25) and (26), so it does not appear explicitly.
Continuing, we find the k th iterated solution to Eqs. (27) and (28) for F (s, U ) and G(s, U ) by making the substitutions in the right-hand sides of Eqs. (31) and (32) and replacing F 1 (s, U ) and G 1 (s, U ) on the left-hand side by F k+1 (s, U ) and G k+1 (s, U ). The resulting expressions for F k+1 and G k+1 add new terms proportional to a k 1 , which again are rational functions of U , with denominator of higher power in U than the numerator. All the terms in F k+1 (s, U ), G k+1 (s, U ) decrease at least as rapidly as 1/|U | for |U | → ∞, and the only singularities are poles at known locations, so we can again calculate the Laplace inversion from U space to τ space analytically. At each stage, we can write the inverse transforms as with the functions k(a 1 , b 1 , s, τ ) expressed as power series in the NLO expansion parameter a 1 whose coefficients are analytic functions of s and τ . These expressions rapidly become too complicated and too lengthy to reproduce here, but are easily calculated using a program such as Mathematica [14]. After numerical Laplace inversion [1,2] of the k's from s to v space, suppressing their explicit dependence on a 1 and b 1 , we define their Laplace inverses as so that we can write the decoupled solutions in (v, τ ) space as the convolutionŝ Finally, recalling that v ≡ ln(1/x), we can transform the above solutions back into the usual space, Bjorken-x and virtuality Q 2 , enabling us to write the NLO decoupled solutions, F s (x, Q 2 ) and G(x, Q 2 ), which require only a knowledge of F s0 (x) and G(x) at Q 2 0 , where evolution is started. In order to insure continuity across the boundaries Q 2 = M 2 c and M 2 b , we first evolve from Q 2 0 (where, e.g., Q 2 0 = 1 GeV 2 for the MSTW group [3]) to M 2 c and use our evolved values ofF s0 (v) andĜ 0 (v) for a new starting values of F s0 (v) andĜ 0 (v). We then evolve to M 2 b , repeating the process, thus insuring continuity of F s (x, Q 2 ) and G s (x, Q 2 ) at the boundaries where n f changes.

III. NON-SINGLET SECTOR
For non-singlet distributions F ns (x, Q 2 ), such as for valence quarks, D val = x d(x, Q 2 ) −d(x, Q 2 ) -the difference between quark distributions-we can schematically write the logarithmic derivative of F ns as the convolution of F ns (x, Q 2 ) with the non-singlet splitting functions, P LO,ns qq (x) and P N LO,ns qq (x), for LO and NLO, respectively, (using the convolution symbol ⊗), i.e., After changing to the variable v = ln(1/x) and the variable τ , we write The comments that we made in Sec. II about integrals that involve the distribution 1/(1 − e −v ) + also apply here. Going to Laplace space s, we obtain a linear differential equation in τ for the transform f ns (s, τ ). This has the simple solution where and We note that in LO, Φ ns (s) = Φ LO f (s), where Φ LO f (s) has been written out explicitly in Eq. (15). Again, the evaluation of Φ N LO ns (s) is straightforward, but too lengthy to be shown here. We can find any non-singlet solution, F ns (x, Q 2 ), by using the non-singlet kernel K ns (v) ≡ L −1 e τ Φns(s) ; v in the Laplace convolution relationF and then returning to (x, Q 2 ) space. In order to insure the continuity of F ns (x, Q 2 ) where n f changes, we renormalize the starting valuesF ns0 (v) at the boundaries M 2 c and M 2 b , as described previously in a similar context for singlet distributions.

IV. CONCLUSIONS
For the singlet sector of pQCD, we have solved the coupled NLO DGLAP equations and found NLO decoupled analytic solutions for F s (x, Q 2 ) and G(x, Q 2 ), an extension of our earlier work for LO [13]. All that is required is knowledge of the initial distributions F s0 (x) and G 0 (x), at Q 2 = Q 2 0 , where Q 2 0 is the starting value for the evolution. For the non-singlet sector, we have successfully solved the NLO evolution equation for F ns (x, Q 2 ), again in terms of F ns0 (x), the value of F ns (x, Q 2 ) at Q 2 0 . We illustrated this numerically for NLO, calculating the non-singlet valence quark distributions U val = x u(x, Q 2 ) −ū(x, Q 2 ) and D val = x d(x, Q 2 ) −d(x, Q 2 ) for a very large range of x and Q 2 , in excellent agreement with the NLO published MSTW [3] values. We note that these techniques can be extended to arbitrary order in the strong coupling constant α s (Q 2 ), for both the singlet and non-singlet sector.
The results presented here are basically analytic, thus eliminating the need for simultaneous numerical solutions of the singlet and non-singlet DGLAP equations on a two-dimensional lattice in x and Q 2 . They provide new tools for studying pQCD; for example, they can be used to examine directly the sensitivity of an individual evolved distribution to the assumed shapes of its starting distribution. In the future, we hope to apply these techniques to a global fit of experimental F γp 2 (x, Q 2 ) data to determine in LO, a gluon starting distribution, and in NLO, approximate F s and gluon starting distributions. Since these starting distributions will be determined by experimental data, they will be free of predetermined shape hypotheses; this will allow us to find new gluon distributions that are critically needed for the interpretation of results from the Large Hadron Collider.

V. ACKNOWLEDGMENTS
The authors would like to thank the Aspen Center for Physics for its hospitality during the time parts of this work were done. P. Ha would like to thank Towson University Fisher College of Science and Mathematics for travel support. D.W.M. received travel support from DOE Grant No. DE-FG02-04AR41308.