A semi-analytical x -space solution for parton evolution — Application to non-singlet and singlet DGLAP equation

: We present a novel semi-analytical method for parton evolution. It is based on constructing a family of analytic functions spanning x -space which is closed under the considered evolution equation. Using these functions as a basis, the original integro-diﬀerential evolution equation transforms into a system of coupled ordinary diﬀerential equations, which can be solved numerically by restriction to a suitably chosen ﬁnite sub-system. The evolved distributions are obtained as analytic functions in x with numerically obtained coeﬃcients, providing insight into the analytic behavior of the evolved parton distributions. As a proof-of-principle, we apply our method to the leading order non-singlet and singlet DGLAP equation. Comparing our results to traditional Mellin-space methods, we ﬁnd good agreement. The method is implemented in the code POMPOM in Mathematica as well as in Python .

The inner structure of hadrons, the bound states of quantum chromo dynamics (QCD), has been under investigation for the past 50 years [1][2][3][4] and continues to be an active area of research to this day [5][6][7].Since the low energy behavior of QCD is non-perturbative, the structure of hadrons cannot be calculated by perturbative methods.Therefore, at present, parton distributions parametrizing the hadron structure have to be extracted from fits to experimental data [8][9][10].In parallel, there is an ongoing effort to calculate aspects of the hadron structure from first principles employing lattice methods [11].This approach is well suited for the prediction of various static properties of hadrons like the magnetic moment [12], the axial charge [13], and form factors [14][15][16][17][18]. Also (low) Mellin-moments of parton distributions are accessible on the lattice, which have the potential to constrain global parton distribution fits [19].For reviews on recent accomplishments regarding the extraction of parton distributions on the lattice see references [20][21][22].
The distribution functions describing the hadron structure enter observables, e.g.structure functions of deep inelastic scattering (DIS), through factorization theorems [23][24][25] who schematically have the form Here σ h (Q, x) is a hadronic observable depending on a hard scale Q, which is much larger than the non-perturbative QCD mass scale Λ QCD , and other kinematic quantities denoted by x, e.g.transverse momenta or rapidities.σ h (Q, x) can be factorized into a partonic hard part H p and a parton-in-hadron distribution function f p/h up to corrections of higher orders in the scale Λ QCD /Q.The hard part only depends on short-distance physics which can be described perturbatively in terms of partons p.The distribution functions f p/h , which only depend on the universal, i.e. process-independent long-distance physics, can be viewed as probability densities of finding a parton p in the hadron h.They are determined by fitting to available data, see [26] for a review on the state of the art method using neural networks.
Both quantities H p and f p/h depend on the factorization scale µ F separating perturbative and non-perturbative regimes.To all orders in perturbation theory in the -also scale-dependent -strong coupling α s (µ 2 ), the overall scale dependence would cancel in σ h .Truncating the series for H p at a fixed order however results in a residual dependence on µ F .From demanding 0 = d dµ F σ h and knowing the scale dependence of H p , one can derive evolution equations for the scale dependence of the distribution functions f p/h , schematically df p/h dµ F = P ⊗ f p/h . (1.2) Note that even though the distribution functions are intrinsically non-perturbative objects, the evolution equations have perturbatively calculable kernels P .The most basic non-perturbative distribution functions describing the hadronic structure are parton distribution functions (PDFs) [27].They parametrize the probability of finding a parton of longitudinal momentum fraction ξ in a parent hadron.Different PDFs exist for unpolarized [27], longitudinally polarized [8], and transversely polarized hadrons [28].Their evolution equation is the so-called Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation [29][30][31][32], which will be discussed in more detail below.
More general parton distributions, whose evolution has been studied in the literature [33,34], are the so called transverse momentum dependent parton distributions (TMDs) [11,35] and generalized parton distributions (GPDs) [36].Both carry information about the three dimensional structure of the hadron.TMDs depend on the transverse momentum of the parton and are necessary to capture the physics in scattering regimes of low transverse momentum while GPDs carry information about the transverse position of partons within the hadron.
In this paper, we only explicitly discuss the DGLAP evolution of PDFs.The motivation behind this study is to investigate a semi-analytical method that may be useful for more involved equations like the ETQS evolution.For PDF evolution, we do not strive to compete with existing well-established and very successful methods, but rather take the DGLAP equation as a proof-of-principle for our new method, where we can easily compare its advantages and drawbacks relative to existing methods along with its numerical performance.
There are quite a few solution methods for the DGLAP equation in the literature.The common approaches can be broadly divided into two classes, Mellin-space methods and xspace methods.Mellin methods employ the convolution structure of the DGLAP equation, which under Mellin transformation becomes a simple multiplication.This allows for an analytic solution in Mellin-space.The back transformation into x-space is usually performed numerically.Publicly available algorithms based on this idea are partonevolution [47], QCD-PEGASUS [48], and EKO [49].The x-space methods usually employ a discretization in x together with a local or global interpolation ansatz to bring the DGLAP equation into the form of an ordinary differential equation.Algorithms using this method include QCDnum [50], APFEL [51], HOPPET [52], and ChiliPDF [53].The latter has recently been generalized to encompass double parton distributions [54].Instead of using a basis of interpolating functions one can also use a brute-force discretization in µ F [55,56].Yet another approach that has been explored in the literature is expanding the PDFs and evolution kernels in Laguerre polynomials in ln(x) [57].This transforms the integro-differential equation to a summation of a finite number of Laguerre coefficients.The evolution of these coefficients can be determined recursively, a numeric implementation was done in [58].
Our method belongs to the class of x-space approaches.However, instead of discretizing in x we choose ln m (x) m! x n as an overcomplete family of spanning functions for x-space which is constructed such that it is closed under DGLAP evolution.This transforms the DGLAP equation into an ordinary differential equation for the coefficients.Upon truncation to a finite subsystem, the differential equation can be solved by an ordinary matrix exponential, where the matrix in the exponential follows from a Magnus expansion [59].The goal is to use spanning functions which capture the analytic behavior in as few functions as possible such that we can restrict the infinite function basis to a finite set for subsequent numerical treatment with only small errors.We find that to achieve a relative accuracy at the 0.01% level for DGLAP evolution the chosen family requires O(200) members.The choice of spanning functions as well as the selection of a finite subset is anything but unique and might be refined in subsequent studies.
The novel idea behind our approach is to get a grasp on the analytic behavior of PDFs in x-space generated by DGLAP evolution.This sets our method apart from previously applied x-space methods which use discretization in x and interpolate the resulting PDFs.
We think that the approach presented in this work might help to get insight on distribution functions whose behavior is less well-known, such as the T F function.There, the starting point would be to find a suitable basis which is closed under ETQS evolution.Another potential use of our method is the calculation of derivatives of parton distributions, which can be relevant for the study of kinematic limits [60].Due to the analytic form in x, differentiation is straightforward without accumulating uncertainties through numerical differentiation.Furthermore, the proposed method decouples the evolution of the coefficient from the transformation between x-space and Mellin-space which only affects the basis.Hence, it can be applied for the construction of Mellin-space expressions with a simple analytic structure.
This paper is organized as follows.Section 2 summarizes basic notions of the DGLAP equation.In section 3 we introduce our semi-analytical x-space method and apply it to the non-singlet and singlet DGLAP equations, before briefly discussing the relation to Mellinspace.In section 4 we study the numerical performance of our method and compare the results to the Les Houches benchmark PDFs [61] as well as PDFs obtained using established Mellin-space methods.Section 5 concludes our paper and gives a brief outlook on ETQS evolution.Appendix A collects the master integrals needed for applying our method to the leading order (LO) DGLAP equations, appendix B gives the differential equation for the non-singlet coefficients explicitly, while appendix C lists the input distributions we use for numerical evolution.This paper comes with Mathematica and Python implementations of the program POMPOM (Parton evolution with matrices: polynomials, logarithms, and more) which employ the presented method.They can be found as ancillary files on the arXiv.Appendix D gives an estimate of the computation time of the code.

DGLAP equation and evolution basis
In this section, we briefly introduce basic notions of the DGLAP equation [29][30][31][32].At LO, f (x) is the probability density to find a parton f carrying longitudinal momentum fraction x of its parent hadron.This simple interpretation does not strictly hold at higher orders, where the parton densities become renormalization-and factorization-scheme dependent.Through factorization, the PDFs acquire a dependence on the scale µ.How the PDFs change when the scale is varied, is described by the DGLAP equation.An intuitive derivation of the DGLAP equation can be found in standard textbooks [62].More precise arguments based on the operator product expansion and the renormalization group [63,64] are given in [65].
For the physical basis of quark, antiquark and gluon distributions (q i , qi , g), the DGLAP equation reads d d ln µ 2 q i µ 2 , x g µ 2 , x = α s µ 2 2π q j ,q j 1 x dξ ξ P q i q j (ξ) P q i g (ξ) P gq j (ξ) P gg (ξ) The evolution kernels (or splitting functions) can be calculated as perturbative series in the QCD coupling constant α s , Charge conjugation invariance yields P q i q j = P qi qj , P q i qj = P qi q j , ( while due to flavor symmetry, we have P q i g = P qg , P gq j = P qg , and also P q i q j = δ ij P V qq + P S qq , P q i qj = δ ij P V q q + P S q q . (2.4) The flavor-diagonal valence splitting function P V qq starts at order α 0 s , while P V q q and the sea contributions P S qq , P S q q start at order α 1 s .The difference P S qq − P S q q is of order α 2 s [66].Instead of the physical basis, we write the DGLAP equation in terms of the non-singlet and singlet quark PDFs, which are defined as Here, n f is the number of active quark flavors participating in the QCD dynamics.Since the physical basis has 2n f quark distributions, we need 2n f − 1 linearly independent nonsinglet distributions.Note that a linear combination of non-singlet distributions is again a non-singlet distribution.The advantage of this evolution basis is that the non-singlet distributions decouple from the gluon distribution under DGLAP evolution, (2.7) Since the quark splitting functions have some flavor dependence at higher orders, the evolution kernel of the non-singlet equation depends on the specific choice of non-singlet distribution.The scale evolution of the singlet and gluon distributions is governed by where here, (2.9) The LO evolution kernels are ( [62], chapter 4 and references therein) ) As usual, the plus-distribution appearing here is defined through its integral against an arbitrary test function f (ξ) regular for ξ → 1, The next-to-leading order (NLO) evolution kernels can also be found in [62].The complete DGLAP kernels are currently known to NNLO [66][67][68], approximate results for the fourthorder (N 3 LO) quark-gluon splitting function P gq were published recently [69].The running coupling α s is known to N 4 LO accuracy [70][71][72], a review on the empirical and theoretical knowledge of α s can be found in [73].
The momentum fractions carried by the individual partons must add up to the full momentum carried by the hadron, which translates to the momentum sum rule (2.15) Quark number conservation implies the flavor sum rule for valence distributions q v,i ≡ q i − qi As these sum rules hold independent of µ [65], they can be used to check the quality of our evolution.
3 Semi-analytical x-space solution for parton evolution

General idea
To solve a general integro-differential equation of the form where P⊗ denotes an integral operator in x-space, we propose the ansatz where the f m are a suitable set of spanning functions, while all µ-dependence is contained in the coefficients a m .These coefficients are unknown and will be determined by transforming the original integro-differential equation into an ordinary differential equation for the coefficients, The evolution matrix P is obtained by evaluating the integral in eq.(3.1) and collecting the result with respect to the spanning functions.Once determined, the evolution matrix can be applied to any choice of initial condition.Initial conditions in the chosen basis can be determined in several ways.One possibility is a fit of the ansatz to the initial condition given in x-space.If an analytic form of the initial conditions is available, the coefficients in the chosen basis can also be determined by a direct expansion in this basis.This is the method we will employ for the numerical study in section 4, details of our expansion can be found in the accompanying code.If the initial conditions are only known numerically or in case it is preferable for reasons of numerical accuracy or cost, one may employ a suitable interpolation akin to the method proposed in [53].
Eq. (3.3) is formally solved by a matrix exponential.In case P(µ) and P(µ ′ ) do not commute, careful treatment is required.Using the Magnus expansion discussed in the upcoming section 3.2, the solution can still be written in terms of an ordinary matrix exponential.Restricting ourselves to a finite basis, the infinite set of differential equations can be approximated by a finite set, and the matrix exponential can be performed numerically.
Thus the main challenge is finding a suitable set of spanning functions.The initial function f (µ 0 , x) should be well approximated by a small finite set of spanning functions.The coefficients a m (µ 0 ) of this expansion must be small enough, as otherwise small truncation errors in the evolution matrix P may cause large errors.Additionally, the set of spanning functions must be closed under the considered integro-differential equation.While this requirement will in principle be satisfied by a power series in x, a pure power series will often be insufficient as it can lead to large truncation errors.For example, if inserting a pure power series into the integro-differential equation generates powers of ln(x), these logarithms should be included among the spanning functions, since they are only poorly approximated by a truncated power series.This might result in the set of spanning functions to be overcomplete, so the decomposition of f will not be unique, which however does not present an issue.
To illustrate our method and demonstrate its numerical feasibility, we will apply it to the leading order non-singlet and singlet DGLAP equation in this work.In principle, the proposed method is applicable for any evolution equation of the structure of eq.(3.1).

Solving the non-singlet and singlet DGLAP equation
Let us first consider the non-singlet DGLAP equation (2.7).To solve it in x-space, we make the following ansatz, The explicit logarithms in our overcomplete family of spanning functions ln m (x) x n m! are well suited for accurately capturing the small x-behavior, the large x-behavior is accounted for by polynomials.At LO, plugging this ansatz into eq.(2.7) yields where the master integrals I 1 and I 2 are listed in appendix A. The higher-order evolution kernels can be treated in a similar manner.For example, the NLO evolution kernels contain ln(1 ± x), Li 2 (x), which lead to the same master integrals I 1 and I 2 when expanded as a series in x before integration.Inserting the master integrals and collecting with respect to our spanning functions, we compare coefficients on both sides to read off a differential equation for the coefficients a M N µ 2 in matrix form.Since we will need to truncate the series expansion of ln(1 − x) when restricting ourselves to a finite subset of spanning functions, the largest relative errors are to be expected near x = 1.
The differential equation for the coefficients of the LO non-singlet DGLAP equation is explicitly given in appendix B. Structurally, it takes the form Here, M N is the row index while ij is the column index of the matrix P.This matrix equation is solved by Note that when taking higher-order corrections to the splitting functions into account P = P (0) + αs 2π P (1) + αs 2π 2 P (2) + . .., where the different P (n) do not necessarily commute.
Using the Magnus expansion [59], eq.(3.6) can still be solved by an ordinary matrix exponential, such that where ) We see that starting at NNLO, that is O(α 3 s ), commutators of different P (n) contribute.Expressions for Ω n≥3 , which contain nested commutators of the P (n) , can be found in [59] and references therein.Contributions due to Ω 3 are at least of O(α 4 s ).The intermediate points µ k are required if µ is outside the domain of convergence of the Magnus expansion [59].They have to be chosen such that omitted terms in the expansion are sufficiently suppressed, in other words such that the counting of powers of α s in the integrand is not spoiled by a large integration domain.The optimal number of slices K will have to be determined on a case-by-case basis balancing between required accuracy and computation time.A detailed discussion of a Magnus expansion approach towards a differential equation of the form of eq.(3.3) can be found in [59], where its numerical performance is compared to standard Runge-Kutta methods.Of course the latter could also be used for solving eq.(3.6).
We note that using the Magnus expansion in the context of the DGLAP equation has also been proposed recently in [74] to determine an analytical solution in Mellin-space beyond leading order.
The evolution of the singlet quark distribution, which was defined in eq.(2.6), is coupled to that of the gluon distribution and governed by eq.(2.8).To solve the singlet DGLAP equation, we make an ansatz analogous to eq. (3.4), Note that we included x −1 in our set of spanning functions to capture the known, more strongly divergent, small-x behavior of the gluon and sea-quark PDFs.Plugging this ansatz into eq.(2.8), we write the appearing integrals in terms of the master integrals I i defined in appendix A. As in the non-singlet case, we insert these master integrals and collect with respect to our spanning functions, which yields a differential equation for the coefficients, At LO, this equation is solved by We can now take PDFs at an initial scale µ 2 0 , parametrize them in the form of eqs.(3.4) and (3.11), and evolve them to any final scale µ 2 .While in principle, the matrices P have infinite dimensions, in practice the spanning functions ln m (x) x n m! of our ansatz yield only small corrections for large m, n (assuming their coefficients do not grow too fast).Therefore, we can cut off the infinite series in our original ansatz at some M max , N max (m), such that the PDFs are parametrized in the form where ν = 0, −1 in the (non-)singlet case.For such a finite dimensional parametrization, the matrix exponentials in eqs.(3.7) and (3.13) can be solved numerically.Note that we have some freedom to choose which finite subset of basis functions we keep.Generally, the numerical performance will increase by adding more functions, however for a fixed basis size choosing a good cut-off, keeping "important" basis functions in favor of others, will greatly enhance the numerical quality of results.

Relation to Mellin-space
For brevity of notation, we only consider non-singlet PDFs in this section, the observations also apply to the singlet and gluon PDFs.The Mellin transform M with Mellin-space variable s of eq.(3.4) is a mn µ 2 (−1) m (n + s) m+1 . (3.15) As the Mellin transform only acts on the spanning functions and not on the coefficients, the coefficients a mn in Mellin-space are identical to those in x-space.Hence, the scale evolution decouples from the switch between x-space and Mellin-space.The spanning functions in Mellin-space, (−1) m (n+s) m+1 , are meromorphic functions of s on C, with poles on the real axis for s = −n.These properties can for example be exploited to obtain a simple analytic form of the Mellin transform of a PDF given numerically in x-space by using our x-space spanning functions to fit the PDF in the form of eq.(3.14), and afterwards replacing the x-space with the Mellin-space spanning functions.As we will see in the next section, the numerical error in the evolution resulting from the truncation in eq.(3.15) can be controlled.

Numerical analysis
The proposed x-space method to solve eqs.(2.7) and (2.8) is exact for an infinite number of basis functions.The key question is whether we can achieve sufficient accuracy for a realistic finite number of basis functions.The goal of this chapter is to establish the numerical validity of the presented method along with its implementation in POMPOM and to determine a finite set of basis functions giving optimal results.The criteria we will use to assess the quality of the evolution are a comparison with existing methods, consistency under back-evolution, the magnitude of sum rule violations from truncation effects, and convergence when increasing the number of basis functions.
As discussed above, we need to select a finite set of basis functions to solve the infinite systems of coupled differential equations (3.6) and (3.12) numerically.As we will see, choosing different cut-offs has a sizeable effect.Hence, we start the numerical investigation of our method by determining a good cut-off M max , N max (m) for the parametrization in eq.(3.14).Here, we use the Les Houches PDFs as a benchmark for the quality of our method [61].
Having established the validity of our x-space method, we use it to evolve a full set of sufficiently realistic PDFs to various scales between µ 2 = 0.26 GeV 2 and µ 2 = 10 8 GeV 2 .Subsequently, we perform back-evolution as a consistency check and compare the x-space evolution to PDFs obtained using Mellin-space methods.Last, we analyze the numerical convergence properties of our x-space method with increasing number of basis functions.

Choice of basis cut-off
To establish the numerical validity of the proposed method we compare with the Les Houches benchmark PDFs [61].We observe that for a fixed total number of basis functions the quality of the agreement strongly depends on the choice of cut-off.For given M max , we parametrize N max (m) as where ⌊. ..⌋ denotes the floor function.The two parameters α and β allow for the inclusion of proportionately more pure polynomial or pure logarithmic terms and for varying the degree to which we want to exclude products.Here, −1/α is the slope of the cut-off for powers of ln(x) greater than powers of x, −β is the slope of the cut-off for powers of x greater than powers of ln(x).Cut-off functions with different α and β are depicted in figure 1.We observe that the functions ln m (x) m! x n rapidly become very small in the entire domain 0 < x < 1 if both m and n increase simultaneously.This behavior can also be seen in figure 1, where we plotted m! x n as a measure for the overall magnitude of the different basis functions.
However, the optimal cut-off has to take into account not only the magnitude of the basis functions but also the magnitude of the coefficients.These are determined from the initial conditions and the respective elements of the evolution matrix.A typical evolution matrix for a particular choice of total number of basis functions and cut-off is depicted in figure 2. ).The matrix is sparse except for the coupling to logarithm-free terms in the uppermost rows.Most of the entries have magnitude smaller than 1.Except for some densely populated "diagonals", the matrix elements are getting small when approaching the cut-offs.
We observe that the matrix is most densely populated for pure polynomial terms coupling to pure polynomials.With increasing powers of logarithms the coefficients coupling to pure polynomials decrease.The rest of the matrix, which couples terms with at least one logarithm to other terms with logarithms is populated quite sparsely.Only some "diagonals" retain non-zero values.The vast majority of matrix elements has values between −1 and 1 depicted light blue respectively light yellow in figure 2. Only on a few diagonals, due to the presence of harmonic numbers in the master integrals, terms grow logarithmically to values of the order of ±10.Away from these diagonals the values decrease towards the cut-off.These properties can also be read off from the differential equation for the non-singlet coefficients a M N in analytic form given in appendix B.
The key observation which allows for restricting the infinite system to a finite matrix is that the coefficients do not grow strongly.Overall, the numerical performance of the finite system may be spoiled by three sources: • omitting "large" basis functions, • omitting matrix elements which couple strongly to the initial condition, • large coefficients in the initial conditions or a very large number of non-negligible coefficients in the initial conditions.
Note that rescaling the basis functions can shuffle around the overall factors between the three.The only way to avoid large coefficients in the initial conditions is choosing suitable spanning functions, the other two determine what finite cut-off is sensible.
Looking at figures 1 and 2, we see that choosing a finite cut-off drops some basis functions which are not per se small and also omits some not inherently small matrix elements.However, we will nevertheless obtain a precise evolution provided that the omitted basis functions do not couple strongly to the initial conditions.This claim is what needs to be checked numerically in the following.

Comparison to Les Houches benchmark tables
We start the investigation of the numerical importance of the basis cut-off by comparing POMPOM against the Les Houches benchmark tables [61].These are commonly used for checking PDF evolution codes [49].In figure 3 we show the results of the comparison between POMPOM with 200 basis functions and the Les Houches benchmark table for an evolution from µ 2 0 = 2 GeV 2 to µ 2 f = 10 4 GeV 2 in a fixed flavor number scheme with n f = 4.The POMPOM cut-off parameters are chosen as (α, β) = (2,11), this choice will be discussed in more detail below.In the figure on the right we observe that for x ≤ 0.5 the relative difference is less than 10 −4 , for all x ≤ 0.3 there is agreement in all significant digits.For the largest values of x the relative deviations increase due to the smallness of the PDFs near near x = 1.As the figure on the right demonstrates, the absolute error for x ≥ 0.5 is below 10 −6 .Hence we see that even though the cut-off on basis functions in the POMPOM method may lead to large relative deviations in the vicinity of x = 1, the absolute error is under control.
To determine the best cut-off, we optimize the parameters α and β to minimize the deviation from the Les Houches benchmark tables [61].Starting from the quadratic cut-off (α, β) = (0, 0) we perform a gradient descent on the lattice of integer values for α and β until we reach a local minimum.This successive optimization of the cut-off is shown in figure 4.
Testing for 50, 100, 150, and 200 basis functions shows that the optimal cut-off parameters α and β depend on the number of functions.Starting from 100 basis functions there is reasonable agreement for valence-like and sea-like distributions.For 200 basis functions most entries agree to all significant digits.The exception to this are the very small sea-like distributions at x = 0.9.We optimized α and β with respect to the average relative deviation of POMPOM and the benchmark tables.This average is strongly dominated by these points.Flavor and momentum sum rules are violated at a level of below 10 −6 .This gives a quantitative measure for the global accuracy of our method.We observe that optimizing the cut-off improves the POMPOM method by some orders of magnitude.
The optimal cut-off may also depend on the initial conditions.In the following, where we consider a different PDF set, for which no benchmark is available, we will use sum rule violations as a self-contained metric to optimize the cut-off.As seen in figure 4, the cut-off which leads to the best agreement with the benchmark also has low sum rule violation.

Evolution of realistic PDF set
From the comparison with the Les Houches benchmark tables we have seen that POMPOM works properly and preserves sum rules to a high level of accuracy.To further showcase the capability of the proposed new method in the context of DGLAP evolution we change the initial conditions and turn to a more realistic model for PDFs.For the numerical evolution in this section, we use the leading order input distributions at µ 2 0 = 0.26 GeV 2 from [75], which are listed in appendix C and take the running coupling at leading order, which is given by eq.(C.7).Note that while these input distributions have an analytic parametrization in x, discrete input distributions can be interpolated with functions of the form (3.14) and subsequently evolved with the POMPOM method.Technical details of such an interpolation are beyond the scope of the present work.The evolution is performed in the variable flavor number scheme, where one evolves up to a flavor threshold, the resulting distribution is then taken as the new input above the threshold.At leading order there are no non-trivial matching conditions, which only come into play starting at NLO [49,51].Here, the output distribution are immediately suitable to be new inputs without the need of additional steps.This is in contrast to an evolution in Mellin-space where evolution space and output space do not coincide and a transformation is required.
The PDFs from [75] provide us with sufficiently realistic initial conditions which are different from the benchmark to demonstrate that our method does not rely on specifics of the initial distributions.We note that due to the change of initial conditions, the optimal cut-off is also slightly different.As a way to fix the optimal cut-off in the absence of a previously known benchmark, we use the average violation of flavor and momentum sum rules.As discussed in the comparison with the Les Houches benchmarks, these sum rules are a sensible measure for the global quality of the evolution.
The results for 200 basis functions and cut-off parameters α = 4, β = 3 are shown in figure 5.When the scale µ 2 increases, the valence PDFs shift to smaller values of x while preserving flavor number.Simultaneously, the gluon and sea quark distributions strongly increase.Hence, at larger scale, a larger fraction of the proton momentum is carried by the gluons and sea quarks, while the total momentum is preserved.

Consistency of back-evolution
To perform a consistency check of our evolution code, we evolve the initial PDFs from the starting scale µ 2 0 = 0.26 GeV 2 to µ 2 = 10 8 GeV 2 and then back to µ 0 .Of course, for an exact evolution, the back-evolved PDF will exactly coincide with the initial PDF.However, any approximate method may result in errors of the order of the approximations applied.In figure 6 we show the difference between the initial condition expressed in the used basis u v,0 (x) and the back-evolved PDF u v,back (x).Compared to the difference due to truncation between the exact initial condition u v,exact (x) and the basis representation u v,0 (x) we can consider the back-evolution to be exact up to numerical inaccuracies.Absolute differences are plotted rather than relative ones to avoid division by very small values near x = 1.
In our case, due to the evolution employing an exponential, the evolution and backevolution operators are exact inverses up to numerical inaccuracies.Hence, truncation effects of the finite cut-off cancel between evolution and back-evolution.Therefore, on the one hand, back-evolution is not suited to evaluate the numerical impact of cut-off effects, on the other hand exact back-evolution is a theoretically desirable feature.

Comparison to Mellin-space evolution
In this section, we compare POMPOM with evolution in Mellin-space.The initial PDFs are again taken from [75].For the Mellin evolution we use a custom Mathematica code based on the method presented in [48].The Mellin evolution code is included in the Mathematica version of POMPOM for comparison purpose.First, we consider rather small numbers of basis functions.The idea is to investigate how many basis functions are required to obtain qualitatively plausible results.This is potentially interesting if one aims for lowering computational cost on the expense of accuracy, which might be desirable when considering the evolution of multi-parton distributions in the There, the basis needs to span multiple dimensions, resulting in a increased number of required functions.Results for the comparison of POMPOM with different numbers of basis functions to Mellin evolution are given in figure 7.
The evolution up to µ 2 = 10 8 GeV 2 is performed in the variable flavor number scheme.The cut-off parameters are chosen as α = 3 and β = 4, which are found to be the optimal numbers for ∼ 75 basis functions.To showcase the qualitative features, absolute results are plotted on a linear scale.A more detailed look at the level of agreement is taken below.
We see that a sufficient number of basis functions is required to obtain qualitatively satisfactory results.As long as basis functions essential for the description of the initial distribution are omitted from the basis, results are far off from the Mellin values.This number is observed to be ∼ 50 for valence-like distributions and ∼ 100 for sea-like distributions.This difference in required basis size is explained by the additional 1 x terms to describe stronger divergence of the sea distributions at small-x.Now we take a closer look at the level of agreement between POMPOM and Mellin evolution.For this purpose we take a moderately higher number of 200 basis functions with cut-off parameters (α, β) = (2, 2).Evolution is performed to µ 2 = 10 8 GeV 2 as above, the relative and absolute differences between POMPOM and Mellin evolution are plotted in figure 8. To visualize the small-x behavior, the left figure shows relative differences on a logarithmic scale.For the non-singlet evolution we find agreement below the ∼ 10 −5 level for 10 −7 x 0.2 with minimal relative difference around x = 10 −4 .For the singlet evolution the agreement is on this level for x 0.1 and is roughly constant down to x = 10 −8 .In the right figure, we look at the large-x behavior on a linear x-scale.To avoid division by the small values of PDFs near x = 1 absolute differences between POMPOM and Mellin Evolution are presented in this case.For large x the absolute differences stay below 10 −5 and do not increase when approaching x = 1 confirming the findings from the comparison with the Les Houches benchmark in figure 3.
These comparisons show that POMPOM is in reasonable agreement with well-established methods for DGLAP evolution at leading order, demonstrating that the cut-off effect can be controlled with a numerically manageable number of basis functions.Since higher orders in α s will presumably, due to their suppression in α s , only result in moderate corrections to the matrix elements used in POMPOM evolution, one can hope for a comparable level of numerical agreement for a certain number of basis functions.A detailed comparison to existing evolution codes such as HOPPET, APFEL, and EKO is beyond the scope of this paper since it will be meaningful only after POMPOM has been upgraded to a state-of-the-art DGLAP evolution code.Here we only aim for a proof-of-concept of the proposed method to justify its future developments for higher order DGLAP evolution or structurally similar evolution equations.

Convergence for increasing number of basis functions
In this section, we take a closer look at the convergence of POMPOM when increasing the number of basis functions.It is clear that increasing the number of basis functions improves the result.However, how quick the convergence towards the true value turns out to be requires investigation.
When comparing with Mellin-space evolution in the previous section, we have already seen some quantitative examples of convergence.Here, we take a closer look at the quark singlet distribution, which is especially difficult to describe due to its strong divergence for small x and particularly its very small values for x ≈ 1.
We fix the cut-off parameters at α = β = 2 and perform an evolution to µ 2 = 10 8 GeV 2 with the PDF set from from [75] while successively increasing the number of basis functions.Then we compare the evolved PDFs for different numbers of basis functions.
In figures 9 and 10, respectively, we compare the absolute and relative differences for evolutions with between 78 and 269 basis functions in steps of ∼ 40.We observe that the differences between steps becomes smaller and smaller indicating convergence towards zero.The absolute differences are largest for small x since the singlet PDF grows strongly in this region.The lowest absolute differences occur at large x ≈ 1.The relative differences are roughly constant for many orders of magnitude of x between 10 −7 and 10 −1 .For large x ≈ 1 the relative difference deteriorates due to the smallness of the singlet PDF in this region.Overall, we observe satisfactory convergence for a modest linear increase in the number of basis functions.By adjusting the cut-off for higher numbers of basis functions, the rate of convergence could be further improved.

Conclusion and Outlook
We have demonstrated that the proposed x-space method implemented in POMPOM is capable of accurately solving the DGLAP equation for the non-singlet and singlet cases.We obtained good agreement with both the Les Houches benchmark set as well as with Mellin-space evolution of the PDFs from [75].Sum rule violations drop to levels below 10 −6 for O(200) basis functions, relative corrections when further increasing the number of basis functions are below 10 −4 except very close to x = 1.POMPOM is accurate over a wide range of x as well as µ, we tested for 2 Λ QCD ≃ µ < 10 TeV and 10 −7 < x < 1.The largest deviations are observed for x ≈ 1 and originate from truncating of the series for ln(1−x).One could try to improve accuracy in the large-x region by amending the basis by ln(1 − x), however since we achieved a, for our purposes, satisfactory level of accuracy, this is beyond the scope of the present paper.The performed numerical studies demonstrate that for the leading order DGLAP equation the cut-off effect of the POMPOM method can indeed be controlled with a numerically manageable number of basis functions.This is an essential foundation for a potential development of POMPOM towards an state-of-the-art DGLAP evolution program that could for example be used for PDF fitting.
Structurally, the DGLAP equations are very similar to the ETQS evolution equation of the twist-3 Qui-Sterman function.In the singlet case the latter is [46] where S ± and F ± are functions of two momentum fractions.The H ± are linear onedimensional integral operators.Expressions for them can be found in the appendix of [44].The non-singlet case decouples similar to DGLAP evolution.The existing evolution code t3evol [46] for eq.(5.1) uses discretization in both µ and x-space.1This code was for example applied in [77].
Having established the validity of the x-space method proposed in this paper for the DGLAP case, generalizing POMPOM to the evolution equation (5.1) is a logical next step.Here, finding a suitable set of basis functions spanning the two-dimensional space of momentum fractions will be key.Along the lines of the presented work, one may start from a polynomial ansatz in both momentum fraction variables and determine which classes of functions are generated by the integration kernels.These functions are subsequently incorporated in the basis if they are not well represented by truncated power series.
From there, selecting an optimal finite subset of basis functions can be expected to be somewhat more challenging than the DGLAP case.Presumably, since there are two momentum fraction variables, more functions will be required for accurate description of the distributions.Hence, it will be more important to find optimal small sets of basis functions to lower the required number.As we have seen in the DGLAP case, for a fixed number of basis functions, the choice of cut-off can lead to order-of-magnitude improvements in accuracy.
An idea worth investigating might be the use of machine learning techniques for the task of cut-off optimization.As a starting point, the impact of a cut-off fitted by a neural net for the DGLAP case could be assessed by comparing with the two-parameter cut-off used in this work.

C Input distributions and parameters
For numerical comparison of our DGLAP evolution with Mellin-space evolution, we use the sufficiently realistic leading order input distributions from [75]

D Computation time
Since the current version of POMPOM is mainly designed for explorative purposes, time performance is not considered the most important aspect of the code.For a potential future upgrade of POMPOM into a code applicable to PDF fits this will be of relevance.One also has to keep in mind that a comparison between evolution programs is only meaningful in reference to a specified task as was dicussed in [49], where the timing of EKO is compared with APFEL.
To get a first estimate on the time required for PDF evolution with POMPOM, we show initialization times and the time for the evolution to a single scale for non-singlet and singlet evolution in figure 11.Initialization has to be performed only once for a fixed number of basis functions with a certain cut-off to calculate the corresponding evolution matrix.The initialization time grows approximately quadratically with the size of the basis.Due to reuse of already calculated matrix elements for lower number of basis functions, the initialization with higher number of basis functions shown in the upper two plots is slightly quicker than it would be in a standalone calculation.The subsequent evolution to a single scale is much faster in comparison and only grows linearly with increasing number of basis functions.For a reasonable number of basis functions, the evolution to a specific scale takes on the order of 1 s on a standard consumer laptop.

1 19 5
Choice of basis cut-off 10 4.2 Comparison to Les Houches benchmark tables 13 4.3 Evolution of realistic PDF set 15 4.4 Consistency of back-evolution 16 4.5 Comparison to Mellin-space evolution 16 4.6 Convergence for increasing number of basis functions Conclusion and Outlook 20 1 Introduction

Figure 1 .
Figure 1.Plot of the contributions of different basis functions to the non-singlet normalization integral.The dashed areas indicate the included basis functions for different cut-offs of the form (4.1) with 150 basis functions each.

Figure 2 .
Figure 2. Magnitude of entries in the non-singlet evolution matrix for (α, β) = (2, 11).The 180 basis functions ln m (x) m! x n are ordered lexicographically by (n, m), i.e. (1, x, x 2 , . . ., ln(x), x ln(x), • • •).The matrix is sparse except for the coupling to logarithm-free terms in the uppermost rows.Most of the entries have magnitude smaller than 1.Except for some densely populated "diagonals", the matrix elements are getting small when approaching the cut-offs.

Figure 3 .
Figure 3.Comparison between POMPOM with 200 basis functions and the Les Houches benchmark table for an evolution from µ 2 0 = 2 GeV 2 to µ 2 f = 10 4 GeV 2 in a fixed flavor number scheme with n f = 4.

Figure 4 .
Figure 4. Plots comparing our method to the Les Houches benchmark tables for different cut-off functions (4.1).The points are labelled with the corresponding cut-off parameter pairs (α, β).The loss plotted in the upper left is the average deviation from the benchmark.Most of the difference stems from the very small sea-PDFs at x = 0.9.The other three plots show the violation of sum rules as a metric for choosing a cut-off which does not rely on a pre-existing benchmark.

Figure 5 .
Figure 5. Evolution of a full set of PDFs in the variable flavor number scheme using POMPOM.For improved visibility, sea-like PDFs are scaled by a factor of 0.1, the gluon PDF is scaled by 0.01.

Figure 6 .
Figure 6.Comparison of back-evolved PDF with initial condition.The red curve shows the difference between the exact initial condition u v,exact (x) and the initial condition expanded in a basis of 200 functions u v,0 (x), which enters the evolution at µ 2 0 = 0.26 GeV 2 .The blue curve shows the difference between u v,0 (x) and the back-evolved PDF u v,back (x), which has been evolved from µ 2 0 = 0.26 GeV 2 to µ 2 = 10 8 GeV 2 and back to µ 2 0 = 0.26 GeV 2 .

Figure 7 .
Figure 7.Comparison of POMPOM with Mellin-space evolution.For increasing numbers of basis functions, POMPOM results converge towards the Mellin results.

Figure 8 .
Figure 8. Relative and absolute differences between POMPOM and Mellin space evolution for nonsinglet and singlet evolution.The left figure shows the relative differences on a logarithmic x-scale to visualize small-x behavior, the right figure shows absolute differences on a linear x-scale to visualize large x behavior.

Figure 11 .
Figure 11.Time required for initialization and evolution of non-singlet and singlet evolution.The timing was performed with Mathematica 13.1 on a ThinkPad T495 with 8 GB of RAM.