ChiliPDF: Chebyshev Interpolation for Parton Distributions

Parton distribution functions (PDFs) are an essential ingredient for theoretical predictions at colliders. Since their exact form is unknown, their handling and delivery for practical applications relies on approximate numerical methods. We discuss the implementation of PDFs based on a global interpolation in terms of Chebyshev polynomials. We demonstrate that this allows for significantly higher numerical accuracy at lower computational cost compared with local interpolation methods such as splines. Whilst the numerical inaccuracy of currently used local methods can become a nontrivial limitation in high-precision applications, in our approach it is negligible for practical purposes. This holds in particular for differentiation and for Mellin convolution with kernels that have end point singularities. We illustrate our approach for these and other important numerical operations, including DGLAP evolution, and find that they are performed accurately and fast. Our results are implemented in the C++ library ChiliPDF.


Introduction
Theoretical predictions at hadron colliders require parton distribution functions (PDFs), which describe the partonic content of the colliding hadrons. PDFs are nonperturbative objects and their exact form is unknown, such that their handling and delivery in practical applications requires approximate numerical methods. Currently available tools use local interpolation over a finite grid of function values, such as splines. For example, the LHAPDF library [1], which has de facto become the standard interface with which PDFs are provided to the community, implements a cubic spline interpolation.
In this paper, we demonstrate a different approach for the numerical representation of PDFs, which is based on a global, high-order interpolation using Chebyshev polynomials, and which allows for significantly higher numerical accuracy at lower computational cost than local interpolation methods. We have implemented this approach in the C++ library ChiliPDF, 1 which is used for most of the numerical demonstrations in the following. We think, however, that the methods presented in our paper are of interest beyond their implementation in a specific software package.
There is a long history of using families of polynomials for handling PDFs or related quantities. Without any claim to completeness, let us mention a few examples. To express PDFs in terms of their Mellin moments, Bernstein polynomials were proposed in ref. [2] and Jacobi polynomials in ref. [3]. An explicit solution of the evolution equations was obtained in refs. [4,5] by expanding PDFs and splitting functions in Laguerre polynomials. More references and discussion can be found in refs. [6][7][8]. To compute analytically the complex Mellin moments of parton luminosities, refs. [9,10] used an expansion in Chebyshev polynomials. The latter also appear in the parameterization of initial conditions in the PDF fits of refs. [11][12][13][14]. We find that all these methods differ quite significantly from each other and from the method to be presented in this work.
Let us emphasize that we are not concerned here with the question of how best to parameterize input PDFs that are fitted to data, nor with the associated systematic error or bias due to the choice of parameterization. For our purposes, we consider the fitted input PDFs to be known "exactly". Our method could be used to address this parameterization issue as well, but we leave this for future exploration.
What we are concerned with is the numerical implementation and handling of PDFs in practical applications, for instance when using PDFs evolved to some scale to obtain quantitative predictions from analytic cross section formulae or with Monte-Carlo event generators. The interpolation of PDFs comes with an inherent inaccuracy that is of purely numerical origin, similar to numerical integration errors. Such inaccuracies should at the very least be small compared to uncertainties due to physics approximations such as the perturbative expansion. But ideally, one would like such inaccuracies to be negligible or at least small enough to be of no practical concern. We will show that this goal can indeed be achieved with Chebyshev interpolation.
Several observations lead us to believe that the performance of available local interpolation methods is becoming insufficient. Generically, the relative accuracy of the interpolation provided by LHAPDF is expected (and found to be) in the ballpark of 10 −3 to 10 −4 . For high-precision predictions this is already close -perhaps uncomfortably closeto the desired percent-level theoretical precision. Moreover, PDFs typically appear in the innermost layer of the numerical evaluation, whose output is then further processed. For example, they enter at the maximally differential level, which is subject to subsequent numerical integrations. Each step comes with some loss of numerical precision, which means the innermost elements should have a higher numerical accuracy than what is desired for the final result. If an integration routine becomes sensitive to numerical PDF inaccuracies, its convergence and hence the computational cost can suffer greatly, because the integrator can get distracted (or even stuck) trying to integrate numerical noise. We have explicitly observed this effect for quadrature-based integrators, which for low-dimensional integrations are far superior to Monte-Carlo integrators, but whose high accuracy and fast convergence strongly depends on the smoothness of the integrand. Another issue is that with increasing perturbative order, the convolution kernels in cross section formulae become more and more steep or strongly localized. This requires a higher numerical accuracy of the PDF, because its convolution with such kernels is sensitive not just to its value but also to its derivatives or its detailed shape. In ref. [15], it was indeed observed that the limited numerical accuracy of PDFs provided by LHAPDF can cause instabilities in the final result.
Many applications require an accurate and fast execution of nontrivial operations on PDFs, such as taking Mellin moments, computing Mellin convolutions with various kernels, taking derivatives, and so forth. A prominent example is DGLAP evolution, which has been extensively studied in the literature, see e.g. refs. [16][17][18][19][20], and for which there is a variety of codes that solve the evolution equations either via Mellin moments [21][22][23] or directly in x space [24][25][26][27][28]. The latter require in particular the repeated Mellin convolution of PDFs with splitting functions, in addition to interpolation in x.
Another important application is the computation of beam functions and similar quantities, which typically appear in resummation formulae. In multiscale problems, beam functions can depend on additional dynamic variables [29][30][31][32][33][34][35][36][37], in which case their evaluation necessitates fast on-the-fly evaluation of Mellin convolutions and in some cases their subsequent DGLAP evolution. A further area is the calculation of subleading power corrections, where the first and second derivatives of PDFs with respect to x are explicitly required [38][39][40][41][42][43][44]. With cubic splines, the first and second derivative respectively correspond to quadratic and linear interpolation and hence suffer from a poor accuracy. Numerical instabilities in derivatives computed from spline interpolants were for instance reported in ref. [45].
Last but not least, for double parton distributions [46] or other multi-dimensional distribution functions, local interpolation methods become increasingly cumbersome due to the large number of variables and even larger number of distributions. By contrast, our global interpolation approach allows for a controllable numerical accuracy with reasonable memory and runtime requirements. This will be discussed in future work.
Clearly, there is always a trade-off between the numerical accuracy of a method and its computational cost, and the accuracy of different methods should be compared at similar computational cost (or vice versa). In our case, the primary cost indicator is the number of points on the interpolation grid, which controls both the number of CPU operations and the memory footprint. 2 As a result of its low polynomial order, local spline interpolation has an accuracy that scales rather poorly with the grid size, so that even a moderate increase in precision requires a substantial increase in computational effort. This can be (partially) offset by improving the performance of the implementation itself, see for example refs. [28,47]. Such improvements can be made for any given method, but they do not change the accuracy scaling of the method itself.
Global interpolation approximates a function over its domain by a single, high-order polynomial. On an equispaced grid, this leads to large oscillations near the edges of the interpolation interval, such that the interpolant never converges and may in fact diverge exponentially. This is well known as Runge's phenomenon [48]. It is often misinterpreted as a problem of polynomial interpolation in general, which may be one reason why local interpolation methods such as splines tend to be preferred. What is perhaps not sufficiently appreciated [49] is the fact that Runge's phenomenon is caused by the use of an equidistant grid and can be completely avoided by using a non-equidistant grid that clusters the grid points toward the edges of the interval. An example for this is Chebyshev interpolation, which leads to a well-convergent approximation, i.e., one that can be made arbitrarily precise by increasing the number of interpolation points. Moreover, approximation with Chebyshev polynomials on a finite interval is closely related to and thus as reliable as the Fourier series approximation for periodic functions. We find that Chebyshev interpolation for PDFs has an excellent accuracy scaling. As we will show, it easily outperforms local spline interpolation by many orders of magnitude in accuracy for the same or even lower number of grid points.
In the next section, we provide some mathematical background on Chebyshev interpolation as we require it. In section 3, we present the global Chebyshev interpolation of PDFs used in ChiliPDF and compare its accuracy with local spline interpolation. We also discuss methods for estimating the numerical accuracy, as well as the basic operations of taking derivatives and integrals of PDFs. In sections 4 and 5, we describe the implementation and accurate evaluation of Mellin convolutions and of DGLAP evolution with ChiliPDF. We conclude in section 6. In an appendix, we discuss the performance of different Runge-Kutta algorithms for solving the evolution equations.

Chebyshev interpolation
In this section, we give a brief account of Chebyshev interpolation, i.e. the interpolation of a function that is discretized on a grid of Chebyshev points. This includes the topics of differentiation, integration, and of estimating the interpolation accuracy. A wealth of further information and mathematical background can be found in ref. [50]. 3 Throughout this section, we consider functions of a variable t restricted to the interval [−1, 1]. The relation between t and the momentum fraction x of a PDF is specified in section 3.
Chebyshev polynomials. The Chebyshev polynomials of the first and second kind, T k (t) and U k (t), are defined by T k (cos θ) = cos k θ , U k (cos θ) = sin(k + 1)θ sin θ (2.1) for integer k ≥ 0. They are related by differentiation as dT k (t)/dt = k U k−1 (t), and they are bounded by |T k (t)| ≤ 1 and |U k (t)| ≤ k + 1. The relations hold for both V = T and V = U . They show that T k (t) and U k (t) are indeed polynomials, which may not be immediately obvious from eq. (2.1). Both families of polynomials form an orthogonal set, i.e. for k, m ≥ 0 they satisfy The first chapters of this book are available on https://people.maths.ox.ac.uk/trefethen/ATAP.

-4 -
where α 0 = 2 and α k = 1 otherwise. For given N , the Chebyshev points are given by They form a descending series from t 0 = 1 to t N = −1 and satisfy the symmetry property t N −j = −t j . The polynomial T N (t) assumes its maxima +1 and minima −1 at the Chebyshev points. We call the set of Chebyshev points a Chebyshev grid. Using eq. (2.1) and expressing sines and cosines as complex exponentials, one readily derives the discrete orthogonality relations where β 0 = β N = 1/2 and β j = 1 otherwise. Notice that the density of Chebyshev points increases from the center toward the end points of the interval [−1, 1]. This feature is crucial to avoid Runge's phenomenon for equispaced interpolation grids, as discussed in the introduction.
Chebyshev interpolation and Chebyshev series. We now consider the approximation of a function f (t) by a finite sum of Chebyshev polynomials T k (t). Using T k (t j ) = T j (t k ) and the first relation in eq. (2.5), one finds that the sum with the interpolation coefficients In other words, p N (t) is the unique polynomial of order N that equals the function f (t) at the N + 1 Chebyshev points t 0 , . . . , t N . We therefore call p N (t) the Chebyshev interpolant.
A sufficiently smooth function f (t) can be expanded in the Chebyshev series whose series coefficients readily follow from the first orthogonality relation in eq. (2.3). Substituting t = cos θ and using eq. (2.1), we recognize in eqs. (2.9) and (2.10) the Fourier cosine series for the function F (θ) = f (cos θ), which is periodic and even in θ. The Chebyshev series on the interval [−1, 1] is thus nothing but the Fourier series of a periodic function in disguise, with the same excellent convergence properties for N → ∞. The Chebyshev series f N (t) is not immediately useful in practice, because its coefficients can only be obtained by explicitly carrying out the integral in eq. (2.10). The Chebyshev interpolant p N (t), however, can be computed very efficiently (see below) and only requires evaluating f (t j ) at the N + 1 Chebyshev points. The key property of interpolating in the Chebyshev points is that the resulting interpolation coefficients c k approach the series coefficients a k in the limit N → ∞. Their precise relation can be found in [50, chapter 4]. It is therefore guaranteed that p N (t) approaches f N (t) and thus f (t) for N → ∞.
Interpolation accuracy. How accurately p N (t) approximates the function f (t) depends on the smoothness of f and its derivatives. We give here a convergence statement that is useful for the interpolation of PDFs. For typical parameterizations, input-scale PDFs are analytic functions of the momentum fraction x for 0 < x < 1, but nonanalytic at x = 1, where they behave like (1 − x) β with noninteger β. We anticipate that this corresponds to a behavior like (1 + t) β at t = −1 when we map an x interval onto a Chebyshev grid in t.
Suppose that on the interval [−1, 1] the function f and its derivatives up to f (ν−1) with ν ≥ 1 are Lipschitz continuous, 4 and that the derivative f (ν) is of bounded variation V . We recall that a function F (t) is Lipschitz continuous on is finite, where in the second step we have split the interval [−1, 1] into K subintervals with boundaries t 0 , . . . , t K+1 such that on each subinterval F (t) has a definite sign. Under these conditions, one has for all N > ν and all t ∈ [−1, 1]. We note that the Chebyshev series has a similar convergence property, which is obtained by replacing p N (t) with f N (t) and 4V with 2V in eq. (2.12).
For the example of a function f (t) = (1 + t) n+δ with integer n ≥ 1 and 0 < δ < 1, the above statement holds with ν = n. The function and its derivatives up to f (n−1) are Lipschitz continuous on [−1, 1]. The nth derivative is proportional to (1 + t) δ and not Lipschitz continuous but of bounded variation. The (n + 1)st derivative is proportional to (1 + t) −1+δ and thus not of bounded variation, because it diverges at t = −1.
It is important to note that eq. (2.12) provides a bound on the maximal absolute interpolation error anywhere in the interval. Often the absolute interpolation error will be much smaller over most of the interval. In the vicinity of a point where f (t) goes to zero, the relative error can still remain large, and the convergence p N (t)/f (t) → 1 for N → ∞ is in general not uniform over the full interval. We will indeed see this for the interpolation of PDFs close to zero crossings or to the end point x = 1.
Barycentric formula. A simple and efficient way to compute the Chebyshev interpolant is given by the barycentric formula where t j denotes again the Chebyshev points, and the barycentric basis functions are given by (2.14) b j (t) is a polynomial of order N , although this is not evident from eq. (2.14). The number of operations for evaluating the barycentric formula scales like N . The formula is found to be numerically stable in the interpolation interval. We note that it is not stable for extrapolating the function f (t) outside this interval [50, chapter 5]. The representation given by eqs. (2.13) and (2.14) is a special case of the barycentric formula for the polynomial L(u) of order N that interpolates a function f (u) given on a set of N + 1 distinct points u 0 , . . . , u N : with basis functions L N is called the Lagrange polynomial for the pairs of values {u j , f (u j )}. We will again use eq. (2.15) below. The simple form of eq. (2.14) comes from the fact that the Chebyshev points yield the very simple weights λ i = β i (−1) i 2 n−1 /n.
Differentiation. Given the Chebyshev interpolant p N (t) for a function f (t), one can approximate the derivative f (t) = df (t)/dt by the derivative p N (t). Note that in general f is not equal to p N at the Chebyshev points. Obviously, one cannot compute the exact values of f (t j ) from the function values f (t j ) on the grid. The derivative p N (t) is a polynomial of degree N − 1 and thus also a polynomial of degree N (with vanishing coefficient of t N ). It is therefore identical to its own Chebyshev interpolant of order N , so we can compute it on the full interval [−1, 1] by the barycentric formula (2.17) To obtain the values of p N (t j ), we take the derivative of eq. (2.6) using T k = k U k−1 . The resulting discrete sums are easily evaluated using eq. (2.1) and expressing the sine function in terms of complex exponentials. We then obtain the relation Note that the matrix multiplication (2.18) maps a vector f (t k ) = const onto the zero vector, as it must be because the derivative of a constant function is zero. Higher derivatives of the Chebyshev interpolant p N (t) can be computed by repeated multiplication of f (t k ) with the differentiation matrix D jk and subsequent application of the barycentric formula. Since each derivative reduces the degree of the interpolating polynomial by 1, the accuracy of approximating f (n) (t) by p one readily obtains the integration rule with weights This is known as Clenshaw-Curtis quadrature. A detailed discussion of its accuracy (and comparison with Gauss quadrature) can be found in ref. [50, chapter 19] and ref. [51].
Interpolation without end points. The interpolant (2.6) is a sum over Chebyshev polynomials T k . Another interpolant can be obtained from the polynomials U k , namely Using the second relation in eq. (2.5) together with sin This means that q N −2 (t) interpolates f (t) on the same Chebyshev points as p N (t), with the exception of the interval end points. Correspondingly, q N −2 has polynomial degree N − 2 rather than N , as is evident from eq. (2.23). We thus have an alternative approximation for f (t), which can be computed from the same discretized function values as those needed for computing p N (t). For sufficiently large N , one may expect that q N −2 approximates f (t) only slightly less well than p N (t), given that its polynomial degree is only smaller by two units. We may hence use |q N −2 (t) − p N (t)| as a conservative estimate for the interpolation error |f (t) − p N (t)|. 5 To evaluate q N −2 (t), we can again use a barycentric formula, namely The corresponding basis functions can easily be obtained from eq. (2.14) by comparing the general expressions (2.15) and (2.16) for the sets of points t j with j = 0, . . . N or j = 1, . . . N − 1. We note that using the barycentric formula (2.26) Formulae for differentiation and integration of f (t) that use the approximant q N −2 are easily derived in analogy to the formulae that use p N . For differentiation, one obtains For integration, one uses for even j, to obtain an open integration rule This is well known as Fejér's second rule, see e.g. ref. [52]. 6 Using eq. (2.32) to estimate the error of Clenshaw-Curtis integration (2.21) is similar to Gauss-Kronrod quadrature [53,54], where for a grid with 2N + 1 points one has an integration rule of order 3N + 1 and a Gauss rule of order 2N − 1, where the latter uses a subset of the points and is used to estimate the integration uncertainty. 7 The difference in polynomial orders between the two rules is hence larger in this case than for the pair of Clenshaw-Curtis and Fejér rules, so that one may expect the error estimate in the latter case to be closer to the actual error. We shall come back to this in section 3.4.
Given the Chebyshev grid (2.4) with N + 1 points for even N , one might also think of estimating the integration or interpolation accuracy by using the subgrid t 0 , t 2 , . . . t N −2 , t N , which has N/2 + 1 points and is again a Chebyshev grid. For the grids and functions we will consider in this work, this would however give a gross overestimate of the actual error, because the interpolant with N/2 + 1 points is significantly worse than the one with N + 1 points. By contrast, using interpolation without the end points of the original N + 1 point grid, we obtain error estimates that are rather reliable, as will be shown in section 3.4.

Interpolation strategy
To interpolate a parton distribution function f (x) in the momentum fraction x, we interpolate the functionf (x) = xf (x) in the variable u = ln x. For a given minimum momentum fraction x 0 , the interval [x 0 , 1] is thus mapped onto the interval [u 0 , 0]. For reasons discussed below, we usually split the x interval into a few subintervals. On each subinterval, we perform a linear transformation from u = ln x to t ∈ [−1, 1] and introduce a Chebyshev grid in t, which is used to interpolate the function as described in section 2. To specify the full grid, which is a conjunction of k subgrids, we use the notation where the x i are the subinterval boundaries and n i = N i + 1 is the number of Chebyshev points for subgrid i. We will refer to this as an (n 1 , n 2 , . . . , n k )-point grid. Note that adjacent subgrids share their end points, so that the total number of grid points is n pts = To be specific, let us consider one subgrid with N + 1 points u 0 , . . . , u N , which is mapped by a linear transform onto the Chebyshev grid t 0 , . . . , t N given by eq. (2.4). The corresponding grid points in x are x i = e u i . We can then interpolate the PDF using the barycentric formulã Here we have used that the form of the barycentric basis functions (2.14) remains unchanged under a linear transform of the interpolation variable. Formulae analogous to eq. (3.2) can be used to interpolate functions derived from PDFs, such as their derivatives (see section 3.3) or Mellin convolutions of PDFs with an integral kernel (see section 4). One reason to use subgrids concerns error propagation. As seen in eq. (3.2), the interpolation off at a certain value x involves the values off on all grid points in the interpolation interval, although the weight of points close to x is higher than the weight of points far away. In an x region wheref is much smaller than its maximum in the interval, numerical errors from regions of largef can strongly affect the interpolation accuracy. This is not much of an issue in our tests below, wheref is computed from an analytic expression, but it does become important when interpolating a PDF that has been evolved to a higher scale and is thus affected by numerical errors from the solution of the DGLAP equations.
So on one hand, the accuracy depends on the behavior of the interpolated function on each subgrid. On the other hand, using multiple subgrids for a fixed total number n pts of points decreases the polynomial degree of the interpolant on each subgrid, and the accuracy quickly degrades when the polynomial degree becomes too small. Hence, for given n pts and x 0 , there is a certain range for the number of subgrids that give the best performance for typical PDFs. We find that the optimum is to take 2 or 3 subgrids for the values of n pts and x 0 used in the following.
To study the accuracy of Chebyshev interpolation for typical PDFs, we consider a number of representative test functions, covering a broad range of shapes and analytic forms, where T k denotes the Chebyshev polynomials defined in eq. (2.1). These functions correspond to PDFs at the input scales of several common PDF sets. Specifically, we have f 1 =ū at NNLO for ABMP16 [55], f 2 = g at LO for MMHT2014 [14], f 3 = g at NLO for HERAPDF2.0 [56], and f 4 = d v at NLO for JR14 [57]. We have studied several other , which decreases very steeply and behaves like a typical gluon density at high scale. The results shown in the following are representative of this more extended set of functions. Throughout this work, the relative numerical accuracy of interpolation for a given quantity is obtained as where the exact result is evaluated using the analytic form of the functions in eq. (3.4). In most plots, we also show the exact result itself as a thick black line, which is solid (dashed) when the result is positive (negative). The relative accuracy for other numerical operations is obtained in full analogy to eq. (3.5).

Interpolation accuracy and comparison with splines
We now compare Chebyshev interpolation with local interpolation. As a prominent example for a method widely used in high-energy physics calculations, we take the interpolation provided by the LHAPDF library [1], which has become a standard interface for accessing parton densities. LHAPDF offers linear and cubic splines in either x or ln x. We use cubic splines in ln x, which is the default in LHAPDF and gives the most accurate results among these options. 8 For brevity, we refer to these as "L-splines" in the following.
Since spline interpolants come in a wide variety, we also consider the cubic splines used by the Interpolation command of Mathematica (versions 11 and 12) and refer to them as "M-splines". For M-splines, both the first and second derivative of the interpolant are continuous, whilst for L-splines only the first derivative is continuous but the second is not.
PDF set x 0 n pts ρ = −n pts / log 10 (x 0 ) The interpolation grids used in LHAPDF depend on the PDF set and cover a wide range in the minimum momentum fraction x 0 and in the total number of grid points n pts , as shown in table 1 for several common PDF sets. The resulting average density of points per decade in x is ρ = −n pts / log 10 (x 0 ) and essentially determines the numerical accuracy of the spline interpolation. For the sake of comparison, we use the grids with the smallest and the largest density among common PDF sets, which happen to be the MMHT2014 grid (n pts = 64 with ρ = 10.7) and the HERAPDF2.0 grid (n pts = 199 with ρ = 33.1).
In figure 1 we compare the spline interpolation on the low-density grid (MMHT2014 grid with n pts = 64) with Chebyshev interpolation on a (32, 32)-point grid, which has nearly the same total number of points (n pts = 63). The M-splines turn out to be more accurate than the L-splines, but this comes at the expense of them being more complex to construct. The Chebyshev interpolation achieves a significantly higher accuracy by several orders of magnitude than either of the splines. This reflects that, contrary to splines, Chebyshev interpolation uses polynomials of a high degree. In figure 2, we compare splines on the high-density grid (HERAPDF2.0 grid with n pts = 199) with Chebyshev interpolation on a (40, 32)-point grid with a total of n pts = 71. Here, even with less than half the number of points, the Chebyshev interpolation achieves several orders of magnitude higher accuracy. This also highlights that the interpolation accuracy scales much better with the number of points for Chebyshev interpolation than for splines.
In figure 3 we compare the accuracy of interpolation on a single Chebyshev grid with the accuracy obtained with the two subgrids used in figure 2. We see that for the same total number of points, interpolation on two subgrids is more accurate, as anticipated above.
We observe in figures 1 to 3 that the relative accuracy varies with x for Chebyshev   interpolation somewhat more than it does for splines. In fact, the absolute accuracy of Chebyshev interpolation varies much less with x, as can be seen by comparing figure 3 with figure 4. The opposite holds for splines, where the relative accuracy shows less variation than the absolute one. This reflects that Chebyshev interpolation is "global" over the full interpolation interval, whilst splines are quite "local" (although the continuity conditions for neighboring splines lead to some correlation over larger distances in x).
As is seen in figures 1 to 3, the relative accuracy degrades in the limit x → 1 for both splines and Chebyshev interpolation. This is not surprising, because in this limit the PDFs in eq. (3.4) approach zero. Moreover, they behave like (1 − x) β with noninteger β and are hence nonanalytic at x = 1. This behavior cannot be accurately reproduced by inter-    polating polynomials in the vicinity of x = 1, whatever their degree. We emphasize that this problem concerns the relative interpolation accuracy. As seen in figure 4, the absolute error of interpolation does remain small for x up to 1, as is expected from eq. (2.12). This is sufficient for many practical purposes, including the evaluation of convolution integrals that appear in cross sections or evolution equations. If high relative accuracy is required at large x, one needs to use subgrids with a sufficient number of points tailored to the region of interest.

Differentiation and integration
We now turn to Chebyshev interpolation for derivatives of the functionf (x) = xf (x). We use the barycentric formula (2.17) and its analog for the second derivative, multiplying of course with the Jacobian for the variable transformation from t to x. For comparison, we also take the numerical derivative (f sp (x + h) −f sp (x − h))/2h of the L-spline interpolant f sp (x). We use a variable step size h = 10 −4 x, having verified that the result remains stable when decreasing h even further. The second derivative off sp is evaluated in an analogous way. In figure 5 we consider the quantities and the accuracy of evaluating them using the two methods just described. We see that with Chebyshev interpolants, a significantly higher accuracy is obtained. This is not surprising, since the polynomials approximating the derivatives are of a high degree in that case, whereas with cubic splines, one locally has a quadratic polynomial for the first derivative and a linear approximation for the second derivative. For the low-density grid with n pts = 64, the L-splines in fact give errors around 10% for the first and 100% for the second derivative in parts of the x range. We note that on each Chebyshev grid, the absolute accuracy of the derivatives (not shown here) has a milder variation in x than the relative one, following the pattern we saw in figure 4 for xf (x) itself.
To explore the accuracy of numerical integration, we consider the truncated Mellin

moments
where the lower integration limit is the lowest point x 0 of the grid. In figure 6, we show the dependence of the relative accuracy on the total number of grid points when using 1 to 4 subgrids of equal size. The advantage of using subgrids is clearly seen, especially for high moment index j, which emphasizes the large x region of the integrand. We find that taking 2 or 3 subgrids gives best results for a wide range of j. The disadvantage of using interpolants with lower polynomial degree becomes the dominant limiting factor with 4 or more subgrids.

Estimating the numerical accuracy
Let us now take a closer look at methods to estimate the numerical accuracy of interpolation or integration with Chebyshev polynomials. An obvious option is to re-compute the quantity of interest with an increased number of points. However, the examples in figure 6 show that the accuracy is not a strictly decreasing function of n pts , and we see fluctuations with local minima and maxima as n pts varies by about 10 units. For a reliable estimate, one should therefore take a sufficiently large increase in n pts . Increasing n pts in several steps provides an additional way of ensuring that the estimate is trustworthy. The procedure just described can give sound accuracy estimates and is what we will adopt for assessing the accuracy of DGLAP evolution in section 5.3. However, since the result must be evaluated multiple times on increasingly dense grids, this procedure can be computationally expensive, especially when the cost scales more than linearly with the number of grid points. This is indeed the case of DGLAP evolution, where the evaluation of convolution integrals involves n pts × n pts matrices.
An alternative with much less computational overhead is to estimate the accuracy from the difference between Chebyshev interpolation and interpolation on the same grid without the end points, as described in section 2. This does not require any additional function evaluations. We compare this estimate with the actual interpolation accuracy in figure 7. In the subinterval for large x, the estimate is very close to the true error, whilst in the subinterval for small to intermediate x it somewhat overestimates the actual numerical error. The amount of overestimation is highest close to the interval limits, which is not surprising, because this is where the two interpolation methods differ most. The corresponding comparison for differentiation is shown in figure 8 and follows the same pattern. We also note that the quality of the estimate for other sample PDFs is as good as or even better than in the examples shown here.
The analog for integration of the previous method is to estimate the accuracy of Clenshaw-Curtis integration from the difference to the Fejér quadrature rule. 9 In figure 9, we compare this estimate with the actual numerical error for the Mellin moment (3.7). We find that it tends to overestimate the actual error by two to three orders of magnitude, which is more than what we saw for interpolation. For comparison, we also  In the accuracy estimates, the results for a given integration rule are first added for the two subintervals, and then the absolute difference between the involved higher-order and lower-order rules is taken.
show the relative accuracy obtained with the widely used Gauss-Kronrod rule, using the same two subintervals and the same number of grid points per subinterval as for Clenshaw-Curtis. As is commonly done, the accuracy estimate for Gauss-Kronrod is obtained from the difference between the nominal Gauss-Kronrod rule and the lower-order Gauss rule on a subset of the grid points. We see that the Gauss-Kronrod rule has the highest accuracy. In practical applications, the true result is however not known, and an integrator is only as good as its accuracy estimation. Here, the accuracy estimate for Gauss-Kronrod is much worse than our estimate for Clenshaw-Curtis. 10 As discussed at the end of section 2, this is not unexpected given the comparison of polynomials orders, but it is interesting to see that the quantitative effect can be as pronounced as in our examples. We also note that the quantitative differences between the integration methods vary significantly with the shape of the integrands, as is evident from the two panels in figure 9.
Overall, as far as numerical accuracy estimates are concerned, using interpolation, differentiation, or integration without the end points of a Chebyshev grid can be considered as an inexpensive estimate of reasonable quality. It is not an alternative to the high-quality estimate one can achieve by a stepwise increase in the number of grid points. Its biggest advantage is thus to provide a fast and reliable indicator whether the accuracy of a result is sufficient or calls for a dedicated, more expensive investigation.

Mellin convolution
In this section, we consider the Mellin convolution of a PDF with an integral kernel, such as a DGLAP splitting function, a beam-function matching kernel, or a hard-scattering coefficient. In terms of the scaled PDFf (x) = xf (x), we wish to compute Here, K(z) is the scaled kernel, given e.g. by K(z) = z P (z) for a DGLAP splitting function P (z). 11 Let us see how eq. (4.1) can be evaluated in a discretized setting. For simplicity, we first consider a single Chebyshev grid for the range x 0 ≤ x ≤ 1 as described in section 3.1. The convolution (K ⊗f )(x) is a function of x over the same domain asf (x). Hence, we can interpolate it in complete analogy tof (x) itself. For this purpose, we only need to evaluate eq. (4.1) at the grid points x i , where in the last step we used the barycentric formula (3.2) to interpolatef (x i /z) under the integral. Introducing the kernel matrix we can compute the function values of K ⊗f on the grid by a simple matrix multiplication, Importantly, the matrix K ij only depends on the given kernel and grid but not onf . It can thus be pre-computed once using standard numerical integration routines.
In general, K(z) is a distribution rather than an ordinary function. Restricting ourselves to plus and delta distributions, we write where K sing (z) is the singular part of the kernel and contains plus distributions. The function K reg (z) contains at most an integrable singularity at z = 1, for instance powers of ln(1 − z). The separation between K sing and K reg is not unique and may be adjusted as is convenient. The convolution (4.1) can now be written as where can be evaluated analytically if the separation between K sing and K reg is chosen appropriately. Due to the explicit subtraction off (x) in the first integral in eq. (4.6), the plus prescription in K sing can now be omitted, because for the logarithmic plus distribution of degree n. However, a PDF vanishes much faster for x → 1 than K d (x) diverges, such that the product K d (x)f (x) tends to zero in that limit. Applying eq. (4.6) to the evaluation of the K ij matrix in eq. (4.3), we have where u i and u j are grid points in u, and where we have changed the integration variable to v = ln z. The plus prescription in K sing can again be omitted, because the term in square brackets vanishes at least linearly in v for v → 0. Hence, all integrals can be evaluated numerically, for which we use an adaptive Gauss-Kronrod routine.
Given that x N = 1, the matrix element K NN is ill defined if K contains a plus distribution, because K d (x) diverges for x → 1 in that case. This does not present any problem: we already noticed that K d (x)f (x) → 0 for x → 1, which translates to K NNfN = 0 in the discretized version. This term may hence be omitted in the matrix multiplication (4.4), along with all other terms K iNfN .
It is not difficult to generalize the preceding discussion to the case of several subgrids in u, each of which is mapped onto a Chebyshev grid. One then has a distinct set of barycentric basis functions for each subgrid. If u j and u i − v are not on the same subgrid, then the basis function b j (u i − v) in eq. (4.10) must be set to zero. As a consequence, the matrix K ij has blocks in which all elements are zero. Let us, however, note that K ij always has nonzero elements below the diagonal i = j and is not an upper triangular matrix.
We now study the numerical accuracy of this method relative to the exact result, which we obtain by the direct numerical evaluation of the convolution integral (4.6). For comparison, we also show the accuracy of the result obtained by approximatingf (x) in eq. (4.6) with its L-splines interpolant and performing the convolution integral numerically. All explicit numerical integrals in this study are performed at sufficiently high precision to not influence the results.
In figure 10, we show the results of this exercise for the convolution of our sample PDFs f i (x) with kernels that have different singular behavior at the end point z = 1. The kernel in the top row is the leading-order DGLAP splitting function P gg , which contains a term L 0 (1 − z). The kernel in the second row is ln 4 (1 − z), which is the most singular term in the three-loop splitting functions P gq and P qg [63,64]. In the bottom row, we take L 5 (1 − z) as kernel, which appears in the N 3 LO corrections to the rapidity spectrum of inclusive Higgs production in pp collisions [15]. In all cases, the accuracy of our Chebyshev based method is several orders of magnitude higher than the one obtained when interpolating PDFs with L-splines and then performing the convolution integral.
The numerical inaccuracies for the convolution with L 5 (1 − z) using L-splines are in the percent range over a wide range of x and much higher than those for the other kernels we looked at. This is in line with the results of ref. [15], which finds that the convolution of L 5 (1 − z) with PDFs from the LHAPDF library can lead to considerable numerical instabilities. The PDFs used in that work are those of the NNPDF3.0 analysis [65], and for the sake of comparison we use the corresponding grid for spline interpolation in the plots on the right of figure 10.
An alternative method to evaluate the convolution is to use the Chebyshev interpolant off and perform the convolution integral itself numerically, analogous to what we did with L-splines above. This avoids the additional interpolation of (K ⊗f ) that happens in the above kernel matrix method. The numerical accuracy of this alternative method differs slightly from the accuracy of the kernel matrix method, but it is of the same order of magnitude. Hence, the additional interpolation of (K ⊗f ) does not introduce a significant penalty in accuracy beyond that of interpolatingf itself. This makes the kernel matrix method far more attractive, as it avoids the evaluation of a numerical integral for each desiredf and x.

Numerical solution of DGLAP equations
In this section, we present our approach to the numerical solution of the DGLAP evolution equations [66][67][68]. Up to order α n+1 s , they read where P (m) is the splitting function at order m. For notational simplicity, we suppress the indices for the parton type and the associated sum on the right-hand side. To solve eq. (5.1), we discretize the scaled PDFsf (x) = xf (x) on a Chebyshev grid and evaluate the Mellin convolutions on the right-hand side as discussed in section 4. The integro-differential equation (5.1) then turns into a coupled system of ordinary differential equations (ODEs), is the kernel matrix for K(z) ≡ zP (m) (z) as defined in eq. (4.3). This system of ODEs can be solved numerically using the standard Runge-Kutta algorithm, which is described in more detail in appendix A. This formulation as a linear system of ODEs by discretization in x is common to many approaches that solve the DGLAP equations in x space [17][18][19][25][26][27]. By using the Chebyshev grid for the discretization in x, the resulting evolved PDF can then be Chebyshev interpolated.
The Runge-Kutta algorithm uses a discretization of the evolution variable t (not to be confused with the argument of Chebyshev polynomials in section 2). It is advantageous if the function multiplyingf on the right-hand side of the evolution equation depends only weakly on t, because this tends to give a uniform numerical accuracy of evolution with a fixed step size in t. We therefore evolve in the variable instead of ln µ 2 . For the running coupling at order n, we write and use a Runge-Kutta routine (see the end of appendix A) to solve for α −1 s as a function of ln µ 2 . At leading order, the DGLAP equation (5.1) then becomes with no explicit t dependence on the right-hand side. Using α s = e −t , we have at NLO where an explicit t dependence appears in the two-loop terms with β 1 and P (1) . The corresponding equations at NNLO and higher are easily written down. The pattern of mixing between quarks, antiquarks, and gluons under DGLAP evolution is well known. Its structure at NNLO is given in refs. [63,69] and remains the same at higher orders. To reduce mixing to a minimum, we work in the basis formed by where q ± = q ±q. Contrary to the often-used flavor nonsinglet combinations u + d − 2s, u + d + s − 3c, etc., the differences between consecutive flavors in eq. (5.7) are less prone to a loss of numerical accuracy due to rounding effects in regions where the strange and heavy-quark distributions are much smaller than their counterparts for u and d quarks.
Note that at NNLO, the combination Σ − has a different evolution kernel than the flavor differences in the second line of eq. (5.7). This is due to graphs that have a t channel cut involving only gluons. We implement the unpolarized splitting functions at LO and NLO in their exact analytic forms. For the unpolarized NNLO splitting functions, we use the approximate expressions given in refs. [63,69], which are constructed from a functional basis containing only the distributions L 0 (1 − z), δ(1 − z) and polynomials in z, 1/z, ln z, and ln(1 − z). With these parameterizations, the numerical evaluation of the kernel matrices in all channels takes about as long for P (2) as it does for P (1) . 12 Both α s and the PDFs depend on the number n F of quark flavors that are treated as light and included in the MS renormalization of these quantities. For the conversion between α s for different n F , we use the matching conditions in ref. [70] at the appropriate order. For the transition between PDFs with different n F , we implement the matching kernels given in ref. [71], which go up to order α 2 s . We have verified that they agree with the independent calculation in ref. [72]. The Mellin convolutions of matching kernels with PDFs are also evaluated as described in section 4.

Validation
To validate our DGLAP evolution algorithm, we compare our results with the benchmark tables in section 1.33 of ref. [73] and section 4.4 of ref. [74]. These tables contain PDFs evolved to the scale µ = 100 GeV from initial conditions given in analytic form at µ 0 = √ 2 GeV. Evolution is performed at LO, NLO and NNLO, either at fixed n F = 4 or with a variable number of flavors n F = 3 . . . 5 and flavor transitions at scales µ c and µ b . The default choice is µ c = m c = √ 2 GeV and µ b = m b = 4.5 GeV. We compare with the LO and NLO results in ref. [73] and with the NNLO results in ref. [74]. The latter are obtained with the same parameterized NNLO splitting functions that we use. Furthermore, the parameterization in eq. (3.5) of ref. [22] is used for the two-loop matching kernel for the transition from a gluon to a heavy quark. We also use this parameterization for the sake of the present comparison, noting that visible differences with the benchmark results appear when we use the exact analytic form for this kernel instead.
The benchmark tables were obtained using two programs, HOPPET [25] and QCD-Pegasus [22], both run with high-precision settings. The former solves the evolution equations in x space, whilst the latter uses the Mellin moment technique. The tables give results for the evolved PDFs at 11 values of x between 10 −7 and 0.9. The results are given with five significant digits, with the exception of several sea-quark combinations at x = 0.9, which are very small and are given to only four significant digits. As specified in ref. [73], evolution with HOPPET was performed by using seventhorder polynomials for the interpolation on multiple x grids spanning the interval [10 −8 , 1] with a total of n pts = 1,170. A uniform grid in ln µ 2 was used, with 220 points between µ = √ 2 GeV and 1000 GeV. The results were verified to be stable within a 10 −5 relative error for x < 0.9 by comparing them with the results of the same program with half the number of points on both the x and the µ grids. Furthermore, they were compared with the results obtained with QCD-Pegasus.
For the comparison with our approach, we use a Chebyshev grid with 3 subgrids [10 −8 , 10 −3 , 0.5, 1] (24,24,24) which has n pts = 70. The DGLAP equations are solved using a Runge-Kutta algorithm with the DOPRI8 method (see appendix A) and a maximum step size in t of h = 0.1. This amounts to about 12 Runge-Kutta steps from the starting scale to µ = 100 GeV. Our results for the benchmark numbers (rounded as stated above) remain the same if we take 40 instead of 24 points for each subgrid in x. They also remain unchanged if we use the same Runge-Kutta method with maximum step size h = 0.3 or h = 0.02.
We compare our results with the numbers reported in tables 2, 3, 4 of ref. [73] and in tables 14 and 15 of ref. [74], which cover unpolarized evolution both with fixed n F = 4 and with variable n F = 3 . . . 5. Whilst the tables also give results for evolution with different scales µ r in α s and µ f in the PDFs, we always set µ r = µ f . We agree with all benchmark numbers, 13 except for those given in our tables 2 and 3. In all cases where we differ, the differences can be attributed to the benchmark results and fall into two categories: 1. The benchmark tables contain a number of entries, marked by an asterisk, for which the results of the two used codes differ in the sixth digit and give different numbers when rounded to five digits. The number with the smaller modulus is then given in the tables. In several cases, our result agrees with number with the larger modulus and in this sense agree with the benchmark results.
2. We differ from the benchmark numbers in five more cases. For the two cases at NNLO, the benchmark results have typographical errors in the exponent. In the other three cases, we differ by one unit in the fifth digit. We contacted the authors of the benchmark tables, who confirmed that indeed their respective codes agree with our numbers [75].  Table 2. Results for evolution with n F = 4 for which we differ from the numbers in the benchmark tables in ref. [73]. The notation for PDF combinations is L ± =d ±ū and q + = q +q for q = s, c.
The number format a b is shorthand for a × 10 b . As discussed in the text, all differences can be attributed to the benchmark results, where entries with an asterisk correspond to category 1 and the entry in blue corresponds to category 2. We note that in the  Table 3. As table 2, but for evolution and flavor matching from n F = 3 . . . 5. The corresponding benchmark tables are in ref. [73] for LO and NLO and in ref. [74] for NNLO.

Numerical accuracy
We now study the numerical accuracy of our method in some detail. We use x 0 = 10 −7 and 3 subgrids. The evolution equations are solved with the DOPRI8 method (see appendix A).
To assess the accuracy, we compare three settings with different numbers of grid points and different Runge-Kutta steps h: To estimate the numerical error due to the discretization in x, we take the difference between settings 1 and 2, whereas the error due to the Runge-Kutta algorithm is estimated from the difference between settings 2 and 3. For the combined error, we take the difference between settings 1 and 3. Note that these estimates correspond to the accuracy of setting 1. We use the same initial conditions at µ 0 = √ 2 GeV as in the benchmark comparison described in the previous subsection. Starting at n F = 3, heavy flavors are added at m c = √ 2 GeV, m b = 4.5 GeV, and m t = 175 GeV. In the remainder of this section, we always evolve and match at NNLO.
The combined discretization and Runge-Kutta accuracy for evolution to µ = 100 GeV and µ = 10 TeV is shown in figures 11 and 12, respectively, both for the individual parton flavors g, q,q, and for the valence combinations q − = q −q. The relative accuracy is better than 10 −7 up to x ≤ 0.8, and much better than that for smaller x. The same holds when we evolve to µ = 1.01m c or µ = 1.01m b , where the charm or bottom distributions are very small. The relative accuracy of u v and d v increases towards small x. This reflects the strong decrease of these distributions in the small-x limit, as we already observed and explained in section 3. The combinations s − , c − , b − , and t − show less variation at small x, and so does their relative accuracy. In figure 13 we show examples for the separate errors due to discretization and the Runge-Kutta algorithm. With our settings, the overall numerical accuracy is entirely determined by the discretization in x, whilst the inaccuracy due to the Runge-Kutta algorithm can be neglected. The Runge-Kutta accuracy for u v is in fact determined by the machine precision except for very large x, as signaled by the noisy behavior of the error curve.
It is often found that backward evolution, i.e. evolution from a higher to a lower scale, is numerically unstable. 14 The structure of the DGLAP equations is such that the relative uncertainties (physical or numerical) of PDFs become larger when one evolves from a scale µ 1 to a lower scale µ 0 . This property is shared by other renormalization group equations in QCD, including the one for α s (µ). To which extent it leads to numerically unreliable results is, however, a separate question.
To study backward evolution within our method, we perform the following exercise. We start with the input PDFs of the benchmark comparison, evolve them with n F = 4 from √ to µ 1 (step 3). We thus have two quantities that are sensitive to the accuracy of backward evolution: • the difference between the output of step 2 and the input to step 1, which corresponds to the evolution path µ 0 → µ 1 → µ 0 , • the difference between the output of step 3 and the input to step 2, corresponding to the evolution path µ 1 → µ 0 → µ 1 .
The relative accuracy for the two evolution paths is shown in figure 14 for the DOPRI8 method with maximum step size h = 0.1. In both cases, we find very high accuracy, in some cases near machine precision as signaled by the noisy behavior of the curves. Relative errors tend to be larger for the evolution path µ 0 → µ 1 → µ 0 , but they are all below 10 −8 for x ≤ 0.9 (except in the vicinity of zero crossings). Note that a high accuracy is also obtained for the small combinations s − , c − , and b − , which are induced by evolution at order α 3 s . We have also repeated our exercise with the widely-used RK4 method and h = 0.03, which requires approximately the same number of function calls as DOPRI8 with h = 0.1 (see the appendix for more detail). This yields relative errors that are orders of magnitude larger for both evolution paths, but still below 10 −5 for x ≤ 0.9 and away from zero crossings. This shows the significant benefit of using Runge-Kutta methods with high order for DGLAP evolution. Nevertheless, even with the standard RK4 method, we find no indication for numerical instabilities of backward evolution in our setup.
Finally, we note that a Runge-Kutta method with less than the 13 stages of DOPRI8 can be useful if one needs to perform evolution in many small steps, for instance when computing jet cross sections with µ ∼ p T for a dense grid in the transverse jet momentum p T . Setting h ∼ 0.03, we find the DOPRI6 method with its 8 stages to be well suited for such situations.  Figure 14. Relative accuracy of NNLO DGLAP evolution at fixed n F = 5 from µ 0 → µ 1 → µ 0 (left) and from µ 1 → µ 0 → µ 1 (right). The scales are µ 0 = 2.25 GeV and µ 1 = 1 TeV, and the DOPRI8 method with maximum step size h = 0.1 is used to solve the evolution equations. The same grid in x is used as for the benchmark comparison in section 5.2.

Conclusions
We have shown that a global interpolation using high-order Chebyshev polynomials allows for an efficient and highly accurate numerical representation of PDFs. Compared with local interpolation methods on equispaced grids, such as splines, our method can reach much higher numerical accuracies whilst keeping the computational cost at a comparable or reduced level.
Not only interpolation, but also differentiation and integration of PDFs are numerically accurate and computationally simple in our approach. The same holds for the Mellin convolution of a PDF with an integral kernel, which can be implemented as a simple multiplication with a pre-computed matrix. In particular, high accuracy is retained if the kernel has a strong singularity at the integration end point, such as ln 5 (1 − z)/(1 − z) + or ln 4 (1 − z). Even for such demanding applications, it is possible to achieve a numerical accuracy below 10 −6 with only about 60 to 70 grid points. Hence, with our method, numerical inaccuracies become safely negligible for practical physics calculations involving PDFs. If desired, the inaccuracy caused by the interpolation can be estimated with modest addi-tional computational effort by using interpolation on the same grid without the end points. This yields an estimate that is appropriately conservative, but not as overly conservative as the estimate we obtain with Gauss-Kronrod quadrature in the case of integration.
By combining our Chebyshev-based interpolation with a Runge-Kutta method of high order, we obtain a very accurate implementation of DGLAP evolution. Using 70 grid points in x and evolving from µ = 1.41 GeV to µ = 10 TeV, we find relative errors below 10 −7 for x between 10 −7 and 0.8. We successfully tested our implementation against the benchmark evolution tables in refs. [73,74]. Performing backward evolution with our method, we find no indications for any numerical instabilities.
Let us briefly point out differences between our implementation of Chebyshev interpolation and the use of Chebyshev polynomials in parameterizations of PDFs. Perhaps most striking is that the order of the highest polynomial used is much larger in our case. There are several reasons for this. First, our approach aims at keeping uncertainties due to interpolation small compared to physics uncertainties. By contrast, as argued in [13], a PDF parameterization need not be much more accurate (relative to the unknown true form) than the uncertainty resulting from fitting the PDFs to data. Second, as we have seen, polynomials up to order 20 or 40 can be handled with ease in our method, whereas PDF determinations naturally tend to limit the number of parameters to be fitted to data. Third, the detailed functional forms used in the PDF parameterizations refs. [11][12][13][14] at a fixed scale µ differ from each other and from the form (3.2) we use for interpolation at any scale µ. It may be interesting to further investigate the impact of these differences on the number of polynomials required for a satisfactory description of PDFs.
Our approach is implemented in the C++ library ChiliPDF, which is still under development and which we plan to eventually make public. To give a sense of its current performance, we note that the evaluation of the interpolated PDF takes less than a microsecond. With the settings in section 5.2, NNLO DGLAP evolution from µ = 1.41 GeV to µ = 100 GeV at n F = 4 takes on the order of 10 to 20 milliseconds. 16 If significantly faster access to evolved PDFs is needed, it may be preferable to pre-compute the µ dependence and use interpolation in both x and µ. This can easily be done using the techniques we have described.
To conclude, let us mention additional features that are already available in ChiliPDF but not discussed here: (i) polarized evolution up to NNLO for longitudinal and up to NLO for transverse and linear polarizations, with flavor matching up to NLO for all cases, (ii) combined QCD and QED evolution of PDFs with kernels up to O(α s α) and O(α 2 ), with QED flavor matching up to O(α), and (iii) evolution and flavor matching of double parton distributions F (x 1 , x 2 , y ; µ 1 , µ 2 ). For the latter, accuracy goals are typically lower than for computations involving PDFs, but the possibility to work with a low number of grid points is crucial for keeping memory requirements manageable. This will be described in a future paper. A broad range of RK methods with different orders and numbers of stages is known [77]. Widely used is the RK4 method, also called the "classic RK method" or simply "the RK method", which has four stages and order p = 4. We also investigate several methods with p > 4, one by Cash and Karp [78] and three by Dormand and Prince [79,80]. According to their order, we denote the latter by DOPRI5, DOPRI6, and DOPRI8. The main parameters of these methods are given in table 4. As noted in ref. [80], the order of the DOPRI6 method is 7 rather than 6 if the function F (t, y) in eq. (A.1) has no explicit t dependence. The numerical study described below suggests that the DOPRI8 method is of order 10 if F (t, y) is t independent. We do not attempt to provide a general proof for this conjecture and hence put a question mark next to the corresponding entry in our table.
To quantify the performance of the different methods for DGLAP evolution, we evolve the sample PDF f 1 (x) given in eq. (3.4) with the DGLAP kernel for the flavor non-singlet combination u + − d + , either at LO or at NLO. As can be seen in eqs. (5.5) and (5.6), LO evolution corresponds to a t independent function F (t, y), whereas NLO evolution implies an explicit t dependence in F (t, y). We evolve from t 0 = 0.7 to t f = 1.7. According to eq. (5.3), this respectively gives α s ≈ 0.5 and α s ≈ 0.18, so that the t dependence of F (t, y) for NLO evolution is quite strong at the lower end of the evolution interval.
To monitor the accuracy of evolution, we consider either the evolved PDF at a given x or the truncated Mellin moment (3.7) with a given moment index j. The results are similar in both cases. In figure 15, we show the accuracy for the moment with j = 3 as a function of N fc = sN steps , where s is the number of stages of the RK method and N steps the number of RK steps between t 0 to t f . Hence, N fc is the total number of calls to the function F (t, y) during evolution. 17 The accuracy is determined from the difference between the solution at the selected N fc and the solution at a value of N fc in the region where the result does not significantly change any more.
For each RK method, we can identify a region in which the accuracy follows a power law (N fc ) −p . We determine the corresponding value p by a fit and find it in reasonable agreement with the order of the method. This indicates that in the corresponding region of N fc the error estimate in eq. (A.4) is indeed applicable.
Our investigation shows the significant advantage of RK methods with high order p for solving the DGLAP equations. For the example shown in the figure, the DOPRI8 method yields a relative accuracy below 10 −8 with a single step (corresponding to 13 function calls), and the limits of machine precision are reached with 50 to 100 function calls.
A similar conclusion can be drawn for the computation of α s (µ) from eq. (5.4). To illustrate this, we evolve α s (M Z ) = 0.118 down to µ 0 = 2.25 GeV with n F = 5. We use the DOPRI8 method and impose a maximum step size h max = 0.2. The evolved value α s (µ 0 ) is then inserted into the exact analytic expression of µ(α s ) at the appropriate perturbative order. Both at NLO and at NNLO, we find that the resulting value of µ agrees with µ 0 at the level of 2 × 10 −15 .