Using analytic models to describe effective PDFs
Salvador A. Ochoa-Oregon David F. Rentería-Estrada Roger J. Hernández-Pinto German F. R. Sborlini Pia Zurita

 

 


 

FeynArts
FeynCalc
Package-X
Mathematica



    #1Eq. (<ref>) UV^a
        ^b
        ^c
        ^d
        ^eIPARCOS-UCM-24-023Parton distribution functions (PDFs) play a pivotal role in hadron collider phenomenology. They are non-perturbative quantities extracted from fits to available data, and their scale dependence is dictated by the DGLAP evolution equations. In this article, we discuss machine learning (ML) strategies to efficiently compute PDFs by bypassing the scale evolution. Analytical approximations to the PDFs as a function of x and Q^2, including up to next-to-leading order (NLO) QCD effects, are obtained. The methodology is tested by reproducing the  set and implementing the analytical expressions in benchmarking codes. It is found that the computational cost is reduced while the precision of the simulations stays well under control.Using analytic models to describe effective PDFs
            
        S. A. Ochoa-Oregon^a,
        D. F. Rentería-Estrada^b,
        R. J. Hernández-Pinto^a,
        G. F. R. Sborlini^c,
        P. Zurita^d,e
            April 29, 2024
        ================================================================================================================
        
        
        
        
        § INTRODUCTION
        
        Breaking the precision frontier in particle physics is a challenging task. The tiny discrepancies among experiments and theoretical predictions might hide new phenomena, and this forces theoreticians to refine as much as possible their methodologies to produce accurate simulations of particle collisions. Even if a plethora of powerful methods are available (see e.g. Ref. <cit.>), most of them require enormous computing resources to achieve the intended precision. In this direction, the purpose of this investigation is to reduce the resource consumption of higher-order calculations for perturbative QFT whilst keeping the precision achieved by the most advanced theoretical results. This will allow to provide fast and reliable results, and also to reduce the environmental impact associated to the high-energy research <cit.>. 
        
        Most of the current higher-order cross-section simulations rely on the factorization theorem <cit.>, i.e., splitting the whole calculation into perturbative and non-perturbative sectors. The process-dependent perturbative contribution is given in terms of well-defined analytical functions (or highly-efficient numerical representations of them). The universal, non-perturbative part, encoded in parton distribution functions (PDFs), can not be computed from fist principles: only their evolution with the factorization scale chosen can be calculated via the DGLAP <cit.> evolution equations. The determination of the PDFs is, instead, achieved by performing global fits to existing data. For practical implementations, the outcome of these fits are given as tables in the relevant kinematic variables with a function that interpolates over them. The interpolation requires to perform several evaluations and the routines are written to be as fast as possible to allow for the efficient computation of physical observables. However, some cross-sections might require a very large number of trials (and running time) to achieve relevant precision, due to the nature of the process. Thus, it is worth exploring the possibility of reducing the CPU/GPU time required by providing an alternative form for the PDFs. Therefore, in this work, we aim to determine an analytical functional form for a set of known PDFs, by approximating their (x,Q^2)-behaviour. This is equivalent to identify approximations to closed analytic solutions of DGLAP equations, a highly non-trivial problem in the context of coupled integro-differential equations <cit.>.
        
        To achieve this ambitious objective, we profit from artificial intelligence. By using Machine-Learning (ML) techniques, we identify suitable functional dependencies to accurately describe the PDF sets. Our thesis is that, avoiding the interpolation over the grids, it will drastically reduce the time required to perform the simulations/computations. Furthermore, the PDFs would be then written in terms of a few parameters, reducing the storage required for them. Also, having access to analytic formulae for PDFs will allow to calculate their derivatives, both in x and Q, which are relevant for performing the matching of F_2 in the threshold regions.
        
        Throughout this work we focus on the set of collinear proton PDFs extracted by the  collaboration, <cit.>, though the methodology is in principle applicable to any other PDF set, provided it is determined via a similar procedure. We stress that our study does not replace, in any way, the need of performing global fits. On the contrary, it is only by having the PDF sets already determined that we can search a posteriori for a functional form.
        
        This work is organised as follows. In Sec. <ref> we present our proposed analytical model for a set of collinear, unpolarised proton PDFs, using  as baseline. The methodology used to determine the functional form is explained in detail in Sec. <ref>. Our results and their comparison with the baseline PDFs can be found in Sec. <ref>. We also check, in Sec. <ref>, the time required to run benchmarking codes using both the available PDF grids and our a posteriori analytical PDFs. We conclude by summarising our findings in Sec. <ref>.
        
        
        
        
        § STARTING HYPOTHESIS
        
        Traditionally, proton PDFs [This also applies, with some differences, to other non-perturbative quantities feasible of phenomenological studies in pQCD, such as fragmentation functions, fracture functions, transverse momentum dependent distributions, etc.] are determined by proposing a functional form in Bjorken x  for the parton densities (or a linear combination of them) at some chosen initial scale Q_0^2. Giving values to the parameters involved and evolving the PDFs to the experimental scales Q^2 using the DGLAP equations, one obtains the distributions that, convoluted with the partonic cross-sections, is compared with data. Repeating the procedure until an adequate description is achieved (according to some chosen criteria), the best fit parameters are determined. This is the case of many sets of proton and nuclear PDFs (e.g. <cit.>), and fragmentation functions (e.g.<cit.>). Other possibility that does not use a proposed parametrization but rather relies on Neural Networks is also employed by the NNPDF collaboration (see e.g. Ref. <cit.> and references within). 
        
        Regarding the parametric form of the distributions, they are usually chosen to be an Euler beta function, with some extra flexibility given by a multiplicative polynomial or exponential. Schematically: 
        
            f_i(x,Q_0^2)=N_i x^α_i (1-x)^β_i P(x,c_ij)  ,
        
        where x is, in the Breit frame, the fraction of momentum of the proton carried by the parton, i indicates a parton flavour or a combination of them (selected to do the evolution), and P(x,c_ij) is some function of x with coefficients c_ij. This form is flexible enough that very different shapes can be achieved by varying the parameters. After the best fit is found, tables in x and Q^2 are made available (in modern times through <cit.>), with fast interpolating routines. 
        
        To determine an analytical expression for the collinear proton PDFs, which we test with the  set at NLO QCD accuracy, we propose that the Q^2 dependence of the PDFs is given by an extension of Eq. (<ref>), with each parameter acquiring a Q^2 dependence. In other words, we start from the assumption that
        
            f_i(x,Q^2)   =N_i(Q^2) x^α_i(Q^2)   (1-x)^β_i(Q^2)
               × P(x,c_ij(Q^2))   ,
        
        with P being, e.g., a polynomial in x. The aim of our work is to determine these functions using ML methods, and explore up to what extent the Q^2 dependence can be reproduced with these functional forms.
        
        
        
        
        § METHODOLOGY
        
        In order to obtain the coefficients from Eq. (<ref>), we performed a two-step fitting procedure. First, we fixed a PDF set within the  framework: in our study, it was , obtained by the  group, which uses solely data from the HERA collider. Then, we generated a grid of 5000-10000 random points in {x,Q} with x ∈ [10^-4,0.65] and Q ∈ [Q_ Min,1000], evaluating the corresponding PDF. The range in x was chosen in accordance to the recommendations done in Ref. <cit.>. For the gluon and light quarks, we set Q_ Min=Q_0≈ 1.37 GeV, whilst we used Q_ Min= 1.5 GeV and Q_ Min= 4.5 GeV for the charm and bottom quark, respectively.
        
        Regarding the parametrization, we draw inspiration from Ref. <cit.>, where at the initial scale Q_0 the authors proposed
        
            x   u_v(x)    =    A_u_v  x^B_u_v (1-x)^C_u_v  (1+E_u_v  x^2)  ,
            
             x   d_v(x)    =    A_d_v  x^B_d_v (1-x)^C_d_v  ,
            
             x  U̅(x)    =    A_U̅  x^B_U̅ (1-x)^C_U̅  (1+D_U̅  x^2)   ,
            
             x  D̅(x)    =    A_D̅  x^B_D̅ (1-x)^C_D̅  ,
            
             x s̅(x)    =    f_s  x D̅(x)   ,
        
        for the valence and sea distributions of light quarks (u, d and s), and 
        
            x   g(x)    =    A_g   x^B_g (1-x)^C_g  - A'_g  x^B'_g (1-x)^C'_g  ,
        
        for the gluon. In these formulae, we are using f_s=0.4 at Q_0 and s ≡s̅, as well as
        U̅ = u̅  ,   D̅ = d̅ + s̅  .
        
        When the scales probed are above the production threshold of a heavy-quark, said quark will be radiatively produced and thus have a nonzero distribution (unless intrinsic heavy flavours are considered, which is beyond the scope of the present work). In particular, in  the optimised general mass variable-flavour-number scheme RTOPT was used to treat charm and bottom PDFs <cit.>. Due to their radiative nature, the gluon has a noticeable impact in the determination of sea and heavy-quark distributions.
        
        Motivated by this, we proposed functional forms like Eq. (<ref>) to fit u_v(x,Q^2) and d_v(x,Q^2), whilst we relied on gluon-like sums of Euler beta functions, i.e. Eq. (<ref>), for the remaining distributions. To perform a quantitative study of the goodness of the approximation, we used the so-called integral error, defined according to
        Δ_i(Q^2) = I[f^ ML_i(x,Q^2),Q^2]-I[f^ HERA_i(x,Q^2),Q^2]/I[f^ HERA_i(x,Q^2),Q^2]  ,
        
        with f_i^ ML and f_i^ HERA being the PDF corresponding to the flavour i obtained with our analytic approximation (ML-PDF from here onwards) and the  set, respectively. The integration operator I[f,Q^2] is given by
        
        I[f,Q^2] = ∫_10^-4^1   dx   f(x,Q^2)   .
        
        We would like to highlight that this definition of the error is suitable to control the validity of the sum rules. Also, since the PDFs contribute to the hadronic cross-section through a convolution with the partonic cross-section, we observed that the integral error successfully provides an error estimation for the physical observables. We will return to this point in Sec. <ref>.
        
        
        
        
        § RESULTS
        
        In this section, we explicitly report the expressions for the analytic approximations to all the PDF flavours provided by  at NLO. We fitted the u_v, u, d_v, d̅, g, s, c and b distributions following the functional forms proposed in each of the following sub-sections. For each value of Q, we extracted the coefficients {A_i,B_i, …} performing a fit in x with the  function <cit.>. Then, we performed a second fit to determine the Q-dependence of these coefficients. In this step, we relied on ML techniques that provide suitable functional forms to build a basis. The functions obtained after this two-step fitting procedure, together with the full set of analytic ML-PDFs are publicly, available in an ancillary file <cit.>. For this reason we refrain from presenting in the manuscript a table with the fitted parameters and functional forms.
        
        
        
        
         §.§ Up quark distributions
        
        The optimal fit was found by considering u_v and u distributions. To this end, we proposed 
        
        x   u_v(x,Q^2)    =    A_u_v(Q^2)   x^B_u_v(Q^2)  (1-x)^C_u_v(Q^2)
           ×   [1+D_u_v(Q^2)x+E_u_v(Q^2)x^2 
        
           +    F_u_v(Q^2)x^3 + G_u_v(Q^2)x^4 + H_u_v(Q^2)x^5]   ,
        
         x   u(x,Q^2)    =    A_u(Q^2)   x^B_u(Q^2)  (1-x)^C_u(Q^2)  ,
        
           ×   [1+D_u(Q^2)x+E_u(Q^2)x^2
        
           +    F_u(Q^2)x^3]   ,
        
        as functional forms to adjust, simultaneously, the dependence on x and Q^2 for u_v and u distributions, respectively. The u̅ distribution is then given by
        u̅(x,Q^2) = u(x,Q^2) - u_v(x,Q^2)   .
        
        
        
        In Fig. <ref>, we present a comparison between our analytic ML-PDF approximations to x  u_v(x,Q^2) (upper plot) and x u̅(x,Q^2) (lower plot) against  PDF set. We appreciate a very good agreement for different values of Q; in particular, we show results for 10 GeV (blue), 100 GeV (green) and 500 GeV (orange).
        
        
        
        
        
         §.§ Down quark distributions
        
        In this case, the optimal fit was achieved by considering d_v and d̅ distributions. Therefore, we wrote:
        
        x   d_v(x,Q^2)    =    A_d_v(Q^2)   x^B_d_v(Q^2)  (1-x)^C_d_v(Q^2)
           ×   [1+D_d_v(Q^2)x^2+E_d_v(Q^2)x^4 . 
        
           +   . F_d_v(Q^2)x^6]   ,
        
         x  d̅(x,Q^2)    =    A_d̅(Q^2)   x^B_d̅(Q^2)  (1-x)^C_d̅(Q^2)
           -   Θ(x_C,d̅-x) A'_d̅(Q^2)   x^B'_d̅(Q^2)  (1-x)^C'_d̅(Q^2)
           ×   [1+D'_d̅(Q^2)x^2+E'_d̅(Q^2)x^4]   ,
        
        with x_C,d̅ = 0.04. We deemed this value to be the best choice to describe the PDFs, by repeating the procedure for a range of x_C,d̅ and selecting the optimal. Then, the d-quark distribution can be recovered by computing
        
        d(x,Q^2) = d_v(x,Q^2) + d̅(x,Q^2)   .
        
        
        We show the comparison between our analytic ML-PDF approximations to x  d_v(x,Q^2) (upper plot) and x d̅(x,Q^2) (lower plot) w.r.t.  PDF set in Fig. <ref>. Again, the agreement is excellent for the complete range of Q considered in this analysis. 
        
        
        
        
         §.§ Gluon distribution
        
        To obtain an approximation to the gluon distribution, which is more involved than the up and down quarks, we proposed a functional form inspired by Eq. (<ref>). Concretely, we used
        
        x   g(x,Q^2) = f_1(x,Q^2) - Θ(x_C,g-x)  f_2(x,Q^2)   ,
        
        with f_i, for i=1,2, as given in Eq. (<ref>). In this way, for values of x below the threshold x_C,g, Eq. (<ref>) offers more flexibility to fit the PDF, in a fully analogous way as done in Ref. <cit.>.
        
        
        
        
        Furthermore, we noticed that the Q-dependence in the coefficient vary significantly in the whole Q range considered (i.e. from Q_0 to 1000 GeV). Hence, we split the analysis in four separate regions:
        
        R_1    =   {Q_0 ≤ Q ≤ 2.5   GeV}  ,
        
         R_2    =   {2.5   GeV≤ Q ≤ 5   GeV}  ,
        
         R_3    =   {5   GeV≤ Q ≤ 150   GeV}  ,
        
         R_4    =   {150   GeV≤ Q ≤ 1000   GeV}  ,
        
        which we determined by examining the Q behaviour of the gluon density.
        Consequently, we defined 
         x   g(x,Q^2)    =    A_g(Q^2)   x^B_g(Q^2)  (1-x)^C_g(Q^2)
           -   Θ(x_C,g-x) A'_g(Q^2)   x^B'_g(Q^2)
           ×     (1-x)^C'_g(Q^2)  ,
        
        for Q ∈ R_1 and 
         x   g(x,Q^2)    =    A_g(Q^2)   x^B_g(Q^2)  (1-x)^C_g(Q^2)
           -   Θ(x_C,g-x) A'_g(Q^2)   x^B'_g(Q^2)
           ×    (1-x)^C'_g(Q^2)[1+D'_g(Q^2) x^2]   ,
        
        for Q ∈{R_2,R_3,R_4}, together with x_C,g=0.1. Again, this last value was fixed by an exploratory procedure.
        
        We emphasize that, since the analysis was done independently for each region, the functions {A_g(Q^2),B_g(Q^2),C_g(Q^2)} and {A'_g(Q^2),B'_g(Q^2),C'_g(Q^2),D'_g(Q^2)} in Eq. (<ref>) have different behaviours in R_2, R_3 and R_4.
        
        In Fig. <ref>, we show our analytic ML-PDF approximation to x  g(x,Q^2) w.r.t. the corresponding gluon PDF from . Even if small fluctuations occur at low-x, the overall agreement is very good. Furthermore, the agreement remains when choosing different values of Q, thanks to splitting the analysis into four regions. 
        
        
        
        
         §.§ Strange quark distribution
        
        For the strange quark we took the usual assumption that s̅≡ s, a common choice for PDFs at NLO that is also imposed for charm and bottom PDFs. After unsuccessfully trying several functional forms based on Euler beta functions, we noticed that the behaviors below and above Q = 4.5 GeV were slightly different. Thus, we used this threshold and split the Q analysis into two regions: 
        
        R_1    =   {Q_0 ≤ Q ≤ 4.5   GeV}  ,
        
         R_2    =   {4.5   GeV≤ Q ≤ 1000   GeV}  .
        
        We defined 
         x   s(x,Q^2)    =    A_s(Q^2)   x^B_s(Q^2)  (1-x)^C_s(Q^2)
           -   Θ(x_C,s(R_1)-x) A'_s(Q^2)   x^B'_s(Q^2)
           ×     (1-x)^C'_s(Q^2)  ,
        
        and 
         x   s(x,Q^2)    =    A_s(Q^2)   x^B_s(Q^2)  (1-x)^C_s(Q^2)
           -   Θ(x_C,s(R_2)-x) A'_s(Q^2)   x^B'_s(Q^2)
           ×    (1-x)^C'_s(Q^2)[1+D'_s(Q^2) x^2]   ,
        
        where x_C,s(R_1)=0.25 and x_C,s(R_2)=0.1 in R_1 and R_2, respectively. Once again, these values were determined by exploring a range of values and choosing the bests ones (i.e., the ones that gave a smaller value of the integral error).
        
        
        
        
        
        In Fig. <ref>, we show our analytic ML-PDF approximation to x  s(x,Q^2) (dashed lines) w.r.t. the corresponding strange PDF from  (solid dots). At low-x (i,.e. below 10^-3), our fit slightly undershoots the prediction from . Still, the agreement is very good for x > 10^-3, and, in particular, for different values of Q. 
        
        
        
        
         §.§ Charm quark distribution
        
        As mentioned above, we considered the heavy quarks as generated radiatively by the gluons, and thus fixed their distributions to zero when the scale is below their respective production threshold. Above it, in the case of the charm quark, we split the analysis in two different regions in Q. Explicitly, we considered: 
        
        R_1    =   {1.47   GeV < Q ≤ 3   GeV}  ,
        
         R_2    =   {3   GeV≤ Q ≤ 1000   GeV}  ,
        
        since m_c^ pole=1.47 GeV according to  fit including NLO QCD corrections <cit.>. Thus, we defined 
         x   c(x,Q^2)    =    A_c(Q^2)   x^B_c(Q^2)  (1-x)^C_c(Q^2)
           ×    [1+D_c(Q^2) x^2]
        
           -   Θ(x_C,c(R_1)-x) A'_c(Q^2)   x^B'_c(Q^2)
           ×     (1-x)^C'_c(Q^2)  [1+D'_c(Q^2) x^2]   ,
        
        and 
         x   c(x,Q^2)    =    A_c(Q^2)   x^B_c(Q^2)  (1-x)^C_c(Q^2)
           ×    (1+D_c(Q^2) x+E_c(Q^2)x^2)
        
           -   Θ(x_C,c(R_2)-x) A'_c(Q^2)   x^B'_c(Q^2)
           ×    (1-x)^C'_c(Q^2)(1+D'_c(Q^2) x^2)   ,
        
        where x_C,c(R_1)=0.1 and x_C,c(R_2)=0.05 in R_1 and R_2, respectively.
        
        We present a comparison between our analytic ML-PDF approximation to x  c(x,Q^2) (dashed lines) versus the values provided by  (solid dots) in Fig. <ref>. The agreement is impressive, both in x and Q. 
        
        
        
        
        
         §.§ Bottom quark distribution
        
        Finally, we studied the bottom quark distribution. Similarly to the case of the charm quark, we started the fitting for values of Q above the mass of the bottom, i.e., Q > 4.5 GeV. We noticed that the quality of our fit significantly increased if we divided the Q range into two regions: 
        
        R_1    =   {4.5   GeV≤ Q ≤ 15   GeV}  ,
        
         R_2    =   {15   GeV≤ Q ≤ 1000   GeV}  .
        
        As we did for the other flavours, we then proposed
         x   b(x,Q^2)    =    A_b(Q^2)   x^B_b(Q^2)  (1-x)^C_b(Q^2)
           ×    [1+D_b(Q^2) x^2]
        
           -   Θ(x_C,b(R_1)-x) A'_b(Q^2)   x^B'_b(Q^2)
           ×     (1-x)^C'_b(Q^2)  [1+D'_b(Q^2) x^2]   ,
        
        and 
         x   b(x,Q^2)    =    A_b(Q^2)   x^B_b(Q^2)  (1-x)^C_b(Q^2)
           ×    [1+D_b(Q^2) x+E_b(Q^2)x^2]
        
           -   Θ(x_C,b(R_2)-x) A'_b(Q^2)   x^B'_b(Q^2)
           ×    (1-x)^C'_b(Q^2)[1+D'_b(Q^2) x^2]   ,
        
        where x_C,b(R_1)=0.1 and x_C,b(R_2)=0.05 in R_1 and R_2, respectively.
        
        As for the other flavours, in Fig. <ref>, we show our analytic ML-PDF approximation to x  b(x,Q^2) (dashed lines) w.r.t. the corresponding bottom PDF from  (solid dots). The agreement is very good for different values of Q, from 10 to 500 GeV. Even if some discrepancies arise for x < 5 · 10^-4, we can state that our formulae globally provides a reliable approximation to the bottom PDF.
        
        
        
        
        § EFFICIENCY AND QUALITY BENCHMARKS
        
        In this section we discuss quantitatively the quality of the ML-PDF approximations found in Sec. <ref> and compare the time required to compute some chosen observables. 
        
        In order to estimate the discrepancies between our analytic ML-PDFs and the original  distributions in a phenomenologically-relevant way, we rely on our definition of integral error given in Eq. (<ref>). In Fig. <ref>, we show that the integral error for almost all the distributions (i.e. u_v, d_v, d̅, c, b and g) is below 0.5  % in the Q ∈ [Q_0,1000   GeV] range. For the strange quark PDF, we appreciate a larger deviation at higher values of Q, reaching up to 1.5  % error for Q ≈ 1000 GeV. Similarly, u̅ has roughly 4  % error for Q < 700 GeV, although it reaches up to 7  % for Q ∈ [780  GeV,1000   GeV]. Despite these discrepancies, in what follows we will demonstrate that their impact in the calculation of physical observables is very limited.
        
        
        
        
        
        Regarding the time needed for the use of the ML-PDFs and , we found a significant difference between simply calling the distributions and actually implementing them in a cross-section calculation. The results presented in the rest of this section were obtained using a desktop PC with a 16-core Intel i7-13700 processor. It is important to remark that nothing has been parallelised. 
        
        
        
         §.§ Run time
        
        
        We start by presenting, in Tab. (<ref>), the run-time difference when performing a call to our PDFs and , for sets of points in the (x,Q^2) space. As we can appreciate, for N_points≤ 10^3-4 the gain is quite substantial. This is due to the fact that calling  has an overhead time for loading and reading the grid, which albeit small, dominates the total run time when one evaluates a small number of points. Therefore, if we were interested in evaluating a handful of points we would be wise in choosing our ML-PDFs over . As we move to larger number of points N_points≈ 10^5 the weight of this overhead starts to dilute, and once we pass the mark N_points≈ 10^8, we reach a region where  the run time of both methods seem to be comparable. As a consequence, it appears that for a code requiring Monte-Carlo integration, it would not make any difference to apply one or the other. We further investigate this in the next sub-section.
        
        
        
        
        
         §.§ Impact in physical observables
        
        Here, we discuss the impact of using our analytic approximations to the PDFs in two realistic cross-section calculations. 
        
        We start by considering pion production in unpolarised proton-proton collisions at √(S_ c.m.)=7 TeV in the central rapidity region (|η^π|<0.5). These simulations make use of a sequential code based on Ref. <cit.>, modified to use the needed PDFs and fragmentation functions (FFs). The Monte-Carlo integrator is , with 10^5 points, the FF set used is <cit.>, and all factorization and renormalization scales have been taken to be equal to the transverse momentum of the pion (denoted p_T). In the upper plot of Fig. <ref>, we present the p_T spectrum using our ML-PDFs (black dots) and  (red line), while the lower plot depicts the ratio between these two quantities. The (light-blue) uncertainty band in the lower plot comes from the Hessian set of theoretical uncertainties of ; it is also present (but hardly noticeable) in the upper plot. From the lower panel, we can appreciate that the difference between both simulations is well below 2  %, and mostly stays within the PDF uncertainty bands. This is much smaller than the typical perturbative error of NLO QCD calculations (generally larger than O(20  %)). Regarding the CPU time, it took 10 minutes 28 seconds using  and 10 minutes 3 seconds with our ML-PDFs, running on the same system described before. From Tab. <ref>, and considering the number of iterations used, one would expect a much larger gain. The outcome of this run serves to emphasize that the matter is far more subtle. It is worth highlighting that this reduction in the computational time is only due to our optimized ML-PDF, since all the other ingredients of the calculation remain the same in both simulations. 
        
        
        
        
        In the second example, we turned to a simulation that required a significantly larger amount of trials to reach a stable output. We tested our ML-PDFs with the calculation of p+p→γ+π including up to NLO QCD effects <cit.>, at √(S_ c.m.)=13 TeV. We consider the central rapidity region, |η^π|,|η^γ|≤ 2.5, imposing p_T^γ∈ [30   GeV,1500   GeV]. We used the same FFs as before, and we took the renormalization and factorization scales to be equal to
        μ = p_T+p_T^γ/2  .
        
        In this case, we chose to run 10^8 points both in order to be in what we expect (from Tab. <ref>) to be a very low time-gain region and because the observable requires a large number of points to produce physical results. We show the outcome in Fig. <ref>. As before, the upper plot presents the cross-section, while the lower plot demonstrates the ratio between the two computations. Given the p_T range explored, there are some oscillations of the output associated to statistical fluctuations. In fact, using 10^7 points gave unphysical results for some p_T values. We can then say that this observable would greatly benefit from using a higher number of points. For our purpose, 10^8 points is enough. In the p_T ≤ 10 GeV region, that dominates the cross-section, the differences are within 5% and are feasible of further reduction. The runtime for  was of 461 minutes, almost 7 hours, while the ML-PDF run took 403 minutes, almost a full hour less. We summarise the runtime of the cross-sections in Tab. <ref>. This gain is much larger than in the previous test, and far more than expected from the table. Given the duration of this run, we refrained from including the PDF theoretical uncertainties. 
        
        
        
        
        To conclude, the results of these two realistic calculations clearly supports the potential of our ML-PDFs to reduce the computational cost of the simulations, keeping under control the uncertainties of the approximation.
        
        
        
        
        § CONCLUSIONS
        
        In this work, we provide an analytic approximation to PDFs using ML-techniques to adjust both their x and Q-dependence. Our starting hypothesis was the assumption that the x-dependence could be reproduced by Eulerian-like functions, whilst all the Q-dependence is embodied within the coefficients.
        
        By doing so, we obtain a reliable approximation to , taking into consideration up to NLO QCD corrections. We show that the integral error is under control for an ample range of x and Q values. In fact, our ML-PDFs were tested for Q ∈ [Q_0, 1000   GeV], presenting (for most of the distributions) deviations w.r.t.  below percent level. It is important to highlight that this is comparable with the error of the PDF sets themselves, and far smaller than the uncertainties introduced by truncating the perturbative expansion of observables in QCD. In fact, we find that the error induced by the ML-PDF in p+p →π 
        is O(1  %) whilst for p+p →γ+π is O(5  %)[In this work, we neglect to study the uncertainty in the parameters of the ML-PDFs. Such a study is beyond the scope of the present work, and we leave it to future studies.]. These errors are much smaller than the O(20-50  %) -or even larger- uncertainties arising from scale variations. Furthermore, we found a non-negligible reduction of the runtime: O(4  %) for p+p →π and O(13  %) for p+p →γ + π.
        
        To conclude, we want to emphasize that the strategy explained in this article is fully applicable to any PDF set, and even to FFs, at any perturbative order (NLO, NNLO and beyond), due to the fact that we made no assumptions about the order at which DGLAP evolution is truncated. Our only hypothesis is that the Q-dependence is fully embodied within the coefficients {A,B,…} described in Sec. <ref>. Following this approach, we can by-pass DGLAP evolution and speed up the calculation of PDFs (and, eventually, FFs). This does not only imply a reduction of the CPU time, but also leads to a reduction of the CO_ 2 footprint and points towards more efficient and sustainable practises in high-energy physics (HEP) <cit.>. For this reason, further developments to obtain analytic approximations to PDFs/FFs are expected and this work can be considered as a first proof-of-concept paving the way for future improvements. As final remark, we want to stress that our approach does not replace the PDF extractions from data; in fact, they constitute an essential input and the basis for ML-PDF determinations as the one presented here.
        
        
        
        
        § ACKNOWLEDGEMENTS
        
        We would like to thank L. Cieri, V. Mateu and G. Rodrigo for fruitful comments about the manuscript. This work is supported by the Spanish Government (Agencia Estatal de Investigación MCIN /AEI/10.13039/501100011033) Grants No. PID2020-114473GB-I00, PID2022-141910NB-I00 and Generalitat Valenciana Grants No. PROMETEO/2021/071 and ASFAE/2022/009. The work of D. F. R.-E. is supported by Generalitat Valenciana (CIGRIS/2022/145).  S. A. O.-O and R. J. H.-P. are funded by CONAHCyT through Project No. 320856 (Paradigmas y Controversias de la Ciencia 2022) and Ciencia de Frontera 2021-2042. R. J. H.-P. is also funded by Sistema Nacional de Investigadores from CONAHCyT. P.Z. acknowledges support, at the early stages of this work, from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Research Unit FOR 2926, grant number 430915485. P.Z. is funded by the “Atracción de Talento” Investigador program of the Comunidad de Madrid (Spain) No. 2022-T1/TIC-24024. The work of G.S. was partially supported by EU Horizon 2020 research and innovation program STRONG-2020 project under grant agreement No. 824093 and H2020-MSCA-COFUND USAL4EXCELLENCE-PROOPI-391 project under grant agreement No 101034371.
        
        
        
        
        
        
        
        
        
        
        30
        fxundefined [1]
         ifx#1
        fnum [1]
          #1firstoftwo
         secondoftwo
         
        fx [1]
          #1firstoftwo
         secondoftwo
         
        
        
        
        
        
        noop [0]secondoftwo
        
        ref[1]@startlink#1@href
        href[1]#1@endlink
        anitize@url [0] `
        12 `$12
          `&12 `#12 `1̂2 `_12 `%12
        startlink[1]
        endlink[0]
        
        rl [1]href #1
        
        
        
        
        
        
        
        
        
        
        
        
        @bib@innerbibempty
        
         [Heinrich(2021)]Heinrich:2020ybq
          authorauthorG. Heinrich,  10.1016/j.physrep.2021.03.006journaljournalPhys. Rept. volume922, pages1 (year2021), http://arxiv.org/abs/2009.00516arXiv:2009.00516 [hep-ph]NoStop
         [Banerjee et al.(2023)Banerjeeet al.]Banerjee:2023avd
          authorauthorS. Banerjeeet al., @noop   (year2023),http://arxiv.org/abs/2306.02837arXiv:2306.02837 [physics.soc-ph]NoStop
         [Collins et al.(1989)Collins, Soper, and Sterman]Collins:1989gx
          authorauthorJ. C. Collins, authorD. E. Soper,
           and authorG. F. Sterman, 10.1142/9789814503266_0001journaljournalAdv. Ser. Direct. High Energy Phys. volume5, pages1 (year1989), http://arxiv.org/abs/hep-ph/0409313arXiv:hep-ph/0409313NoStop
         [Catani et al.(2012)Catani,
          de Florian, and Rodrigo]Catani:2011st
          authorauthorS. Catani, authorD. de Florian,
           and authorG. Rodrigo,  10.1007/JHEP07(2012)026journaljournalJHEP volume07, pages026 (year2012), http://arxiv.org/abs/1112.4405arXiv:1112.4405
          [hep-ph]NoStop
         [Cieri et al.(2024)Cieri,
          Dhani, and Rodrigo]Cieri:2024ytf
          authorauthorL. Cieri, authorP. K. Dhani,   and authorG. Rodrigo,@noop   (year2024), http://arxiv.org/abs/2402.14749arXiv:2402.14749 [hep-ph]NoStop
         [Altarelli and Parisi(1977)]Altarelli:1977zs
          authorauthorG. Altarelli and authorG. Parisi,  10.1016/0550-3213(77)90384-4journaljournalNucl. Phys. B volume126, pages298 (year1977)NoStop
         [Simonelli(2024)]Simonelli:2024vyh
          authorauthorA. Simonelli, @noop   (year2024), http://arxiv.org/abs/2401.13663arXiv:2401.13663 [hep-ph]NoStop
         [Abramowicz et al.(2015)Abramowiczet al.]H1:2015ubc
          authorauthorH. Abramowiczet al. (collaborationH1, ZEUS),  10.1140/epjc/s10052-015-3710-4journaljournalEur. Phys. J. C volume75, pages580 (year2015), http://arxiv.org/abs/1506.06042arXiv:1506.06042 [hep-ex]NoStop
         [Note1()]Note1
          noteThis also applies, with some differences, to other
          non-perturbative quantities feasible of phenomenological studies in pQCD,
          such as fragmentation functions, fracture functions, transverse momentum
          dependent distributions, etc.Stop
         [Bailey et al.(2021)Bailey,
          Cridge, Harland-Lang, Martin, and Thorne]Bailey:2020ooq
          authorauthorS. Bailey, authorT. Cridge,
          authorL. A. Harland-Lang,
          authorA. D. Martin,  andauthorR. S. Thorne,  10.1140/epjc/s10052-021-09057-0journaljournalEur. Phys. J. C volume81, pages341 (year2021), http://arxiv.org/abs/2012.04684arXiv:2012.04684 [hep-ph]NoStop
         [Hou et al.(2021)Houet al.]Hou:2019efy
          authorauthorT.-J. Houet al.,  10.1103/PhysRevD.103.014013journaljournalPhys. Rev. D volume103, pages014013 (year2021),http://arxiv.org/abs/1912.10053arXiv:1912.10053 [hep-ph]NoStop
         [Alekhin et al.(2017)Alekhin, Blümlein, Moch, andPlacakyte]Alekhin:2017kpj
          authorauthorS. Alekhin, authorJ. Blümlein,
          authorS. Moch,  and authorR. Placakyte, 
          10.1103/PhysRevD.96.014011journaljournalPhys.
          Rev. D volume96, pages014011
          (year2017), http://arxiv.org/abs/1701.05838arXiv:1701.05838 [hep-ph]NoStop
         [Helenius et al.(2022)Helenius, Walt, and Vogelsang]Helenius:2021tof
          authorauthorI. Helenius, authorM. Walt,   and authorW. Vogelsang,  10.1103/PhysRevD.105.094031journaljournalPhys. Rev. D volume105, pages094031 (year2022), http://arxiv.org/abs/2112.11904arXiv:2112.11904 [hep-ph]NoStop
         [de Florian et al.(2012)de Florian, Sassot, Zurita, andStratmann]deFlorian:2011fp
          authorauthorD. de Florian, authorR. Sassot,
          authorP. Zurita,  and authorM. Stratmann, 
          10.1103/PhysRevD.85.074028journaljournalPhys.
          Rev. D volume85, pages074028
          (year2012), http://arxiv.org/abs/1112.6324arXiv:1112.6324 [hep-ph]NoStop
         [Eskola et al.(2022)Eskola,
          Paakkinen, Paukkunen, and Salgado]Eskola:2021nhw
          authorauthorK. J. Eskola, authorP. Paakkinen,
          authorH. Paukkunen,  andauthorC. A. Salgado,  10.1140/epjc/s10052-022-10359-0journaljournalEur. Phys. J. C volume82, pages413 (year2022), http://arxiv.org/abs/2112.12462arXiv:2112.12462 [hep-ph]NoStop
         [Duwentäster et al.(2021)Duwentäster, Husová, Ježo,
          Klasen, Kovařík, Kusina, Muzakka, Olness, Schienbein, and Yu]Duwentaster:2021ioo
          authorauthorP. Duwentäster, authorL. A. Husová, authorT. Ježo,
          authorM. Klasen, authorK. Kovařík, authorA. Kusina, authorK. F. Muzakka, authorF. I. Olness, authorI. Schienbein,  and authorJ. Y. Yu, 
          10.1103/PhysRevD.104.094005journaljournalPhys. Rev. D volume104, pages094005 (year2021), http://arxiv.org/abs/2105.09873arXiv:2105.09873 [hep-ph]NoStop
         [Borsa et al.(2022)Borsa,
          Sassot, de Florian, Stratmann, and Vogelsang]Borsa:2022vvp
          authorauthorI. Borsa, authorR. Sassot,
          authorD. de Florian, authorM. Stratmann,  and authorW. Vogelsang, 
          10.1103/PhysRevLett.129.012002journaljournalPhys. Rev. Lett. volume129, pages012002 (year2022), http://arxiv.org/abs/2202.05060arXiv:2202.05060 [hep-ph]NoStop
         [Borsa et al.(2024)Borsa,
          Stratmann, de Florian, and Sassot]Borsa:2023zxk
          authorauthorI. Borsa, authorM. Stratmann,
          authorD. de Florian,  andauthorR. Sassot,  10.1103/PhysRevD.109.052004journaljournalPhys. Rev. D volume109, pages052004 (year2024), http://arxiv.org/abs/2311.17768arXiv:2311.17768 [hep-ph]NoStop
         [Ball et al.(2022)Ballet al.]NNPDF:2021njg
          authorauthorR. D. Ballet al. (collaborationNNPDF), 
          10.1140/epjc/s10052-022-10328-7journaljournalEur. Phys. J. C volume82, pages428
          (year2022), http://arxiv.org/abs/2109.02653arXiv:2109.02653 [hep-ph]NoStop
         [Butterworth et al.(2016)Butterworthet al.]Butterworth:2015oua
          authorauthorJ. Butterworthet al., 
          10.1088/0954-3899/43/2/023001journaljournalJ.
          Phys. G volume43, pages023001
          (year2016), http://arxiv.org/abs/1510.03865arXiv:1510.03865 [hep-ph]NoStop
         [Buckley et al.(2015)Buckley, Ferrando, Lloyd, Nordström, Page, Rüfenacht,
          Schönherr, and Watt]Buckley:2014ana
          authorauthorA. Buckley, authorJ. Ferrando,
          authorS. Lloyd, authorK. Nordström, authorB. Page, authorM. Rüfenacht, authorM. Schönherr,  and authorG. Watt,  10.1140/epjc/s10052-015-3318-8journaljournalEur. Phys. J. C volume75, pages132 (year2015),http://arxiv.org/abs/1412.7420arXiv:1412.7420 [hep-ph]NoStop
         [Thorne(2012)]Thorne:2012az
          authorauthorR. S. Thorne,  10.1103/PhysRevD.86.074017journaljournalPhys. Rev. D volume86, pages074017 (year2012), http://arxiv.org/abs/1201.6180arXiv:1201.6180 [hep-ph]NoStop
         [Research(2008)]reference.wolfram_2023_nonlinearmodelfit
          authorauthorW. Research, @noop titleNonlinearModelFit, howpublished<https://reference.wolfram.com/language/ref/NonlinearModelFit.html>
          (year2008)NoStop
         [Ochoa-Oregon et al.(2024)Ochoa-Oregon, Rentería-Estrada, Hernández-Pinto, Sborlini, and Zurita]ZENODOPDF
          authorauthorS. Ochoa-Oregon, authorD. Rentería-Estrada, authorR. Hernández-Pinto, authorG. Sborlini,  and authorP. Zurita,  10.5281/zenodo.10996239titleAncillary files for "Using analytic models to describe effective
          PDFs" (year2024)NoStop
         [Zhang(2015)]Zhang:2015tuh
          authorauthorZ. Zhang (collaborationH1, ZEUS), 
          10.5506/APhysPolBSupp.8.957journaljournalActa
          Phys. Polon. Supp. volume8, pages957 (year2015), http://arxiv.org/abs/1511.05402arXiv:1511.05402 [hep-ex]NoStop
         [Jager et al.(2003)Jager,
          Schafer, Stratmann, and Vogelsang]Jager:2002xm
          authorauthorB. Jager, authorA. Schafer,
          authorM. Stratmann,  andauthorW. Vogelsang,  10.1103/PhysRevD.67.054005journaljournalPhys. Rev. D volume67, pages054005 (year2003), http://arxiv.org/abs/hep-ph/0211007arXiv:hep-ph/0211007NoStop
         [Hernández-Pinto et al.(2016)Hernández-Pinto, Epele, de Florian,
          Sassot, and Stratmann]Hernandez-Pinto:2016cnc
          authorauthorR. J. Hernández-Pinto, authorM. Epele, authorD. de Florian,
          authorR. Sassot,  and authorM. Stratmann, 
          10.1088/1742-6596/761/1/012037journaljournalJ. Phys. Conf. Ser. volume761, pages012037 (year2016), http://arxiv.org/abs/1609.02455arXiv:1609.02455 [hep-ph]NoStop
         [de Florian and Sborlini(2011)]deFlorian:2010vy
          authorauthorD. de Florian and authorG. F. R. Sborlini,  10.1103/PhysRevD.83.074022journaljournalPhys. Rev. D volume83, pages074022 (year2011), http://arxiv.org/abs/1011.0486arXiv:1011.0486 [hep-ph]NoStop
         [Rentería-Estrada et al.(2021)Rentería-Estrada, Hernández-Pinto, and Sborlini]Renteria-Estrada:2021rqp
          authorauthorD. F. Rentería-Estrada, authorR. J.Hernández-Pinto,  and authorG. F. R. Sborlini,  10.3390/sym13060942journaljournalSymmetry volume13, pages942 (year2021), http://arxiv.org/abs/2104.14663arXiv:2104.14663 [hep-ph]NoStop
         [Note2()]Note2
          noteIn this work, we neglect to study the uncertainty in the
          parameters of the ML-PDFs. Such a study is beyond the scope of the present
          work, and we leave it to future studies.Stop
