Gluons and sea quarks in the proton at low scales

We study the evolution of parton distributions down to low scales by considering several of their Mellin moments. For the initial conditions, we use a broad array of current parton density fits. Confirming earlier findings in the literature, we conclude that current determinations of parton distributions are incompatible with the idea that gluon or antiquark densities are generated by purely perturbative radiation as it is encoded in the DGLAP evolution equations.


Introduction
An outstanding question in QCD is the relation between its basic degrees of freedomquarks, antiquarks, and gluons-and the concept of "constituent quarks", which plays a major role in the description of the spectrum and static properties of hadrons (see section 15 of [1] for a brief overview and references). The internal structure of hadrons at fine spatial resolution is described by parton distribution functions (PDFs) and similar quantities such as transverse-momentum-dependent distributions, generalised parton distributions, or double parton distributions. It is natural to ask how the resulting picture of the proton as a system made of many quarks, antiquarks and gluons can be related with the picture of the proton as a bound state of just two up quarks and one down quark. a e-mail: markus.diehl@desy.de b e-mail: pascal.stienemeier@desy.de A simple and physically intuitive idea, put forward long ago, is that at coarse spatial resolution the proton contains only valence quarks, and that antiquarks and gluons (as well as quarks with low momentum fraction) are generated by radiation, i.e. by splitting processes like q → qg and g → qq that can be computed in QCD perturbation theory [2][3][4]. In a technical implementation of this idea, gluon and antiquark distributions are zero at some low renormalisation scale, and the parton splitting processes encoded in the DGLAP evolution equations generate nonzero gluon and antiquark distributions at higher scales. However, it was already found in the 1990s that PDFs constructed along these lines are in conflict with experimental data [5][6][7]. Modifying the original idea, the Dortmund group of Glück, Reya, and their collaborators performed a long series of PDF fits with a low initial renormalisation scale (typically well below 1 GeV), so that perturbative scale evolution plays a major role in the shape of these distributions at higher scales. This was done both in the unpolarised and in the polarised sector; see [8,9] and references therein.
The simple scenario just described has often been used to compute parton distributions from dynamical quark models; see, for instance, [10,11]. More complicated quantities can be obtained in the same manner, such as generalised parton distributions [12] or double parton distributions [13]. The study in [11] concluded that such approaches "tend to give a qualitative description of the data" but are insufficient at the quantitative level.
Several approaches have been pursued to connect models at low momentum scales with PDFs. In a formulation put forward long ago [14] and used later, e.g. in [15,16], the three "constituent quarks" of a proton have an internal structure that involves gluons and antiquarks. A different line of work is based on the idea that virtual meson fluctuations (often referred to as "meson cloud") provide a natural source of antiquarks in the proton even at a low scale [17]. We refer to section 4.3.1 of [18] for more detail and references, and to [19][20][21] for recent applications. The role of gluons in this context is discussed in [22,23]. Yet a different picture emerges in the chiral quark-soliton model, where the proton has a Dirac sea of quarks and antiquarks at the typical scale μ ∼ 600 MeV of the model [24,25]. As described in these papers, gluon distributions in this model are suppressed parametrically compared with quarks and antiquarks; in physical terms, one has a scenario in which the quark and antiquark degrees of freedom of the model have themselves a structure that gives rise to gluon distributions.
The goal of the present paper is to approach the preceding discussion from a different angle. Starting with our present knowledge of PDFs, as encoded in PDF sets fitted to experimental data by different groups, we wish to investigate the possibility that the distribution functions for gluons or antiquarks, or both, become zero (or at least small) when one evolves them backwards to low scales. At a practical level, evolving PDFs from high to low scales is however delicate, because small changes at the starting scale of evolution quickly blow up. We circumvent this problem by evolving backwards several Mellin moments of parton distributions. The evolution of Mellin moments is described by simple differential equations, which can be solved numerically without numerical stability problems. We can then check the hypothesis that a given PDF becomes zero under evolution to a low scale by checking whether several of its moments evolve to zero at one and the same scale μ. To get a sense of whether perturbative evolution can actually be trusted at low scales, we compare the different orders for which PDF sets are available, i.e. leading order (LO), next-to-leading order (NLO), and next-to-next-to-leading order (NNLO).
We note that the scale evolution of PDF moments, including evolution to low scales, has been considered extensively in spin physics, with a focus on the first moment of the helicitydependent PDFs. Indeed, the "proton spin crisis" was triggered by the observation that the values of these moments extracted from experiment are not consistent with the idea that at a low scale the spin 1/2 of the proton simply arises from the helicities of its three constituent quarks.
For details, we refer to the review [26], and for more recent work on the evolution of spindependent PDF moments to [27][28][29]. In fact, our present work is somewhat similar in spirit to the study in [28]. That work was concerned with the evolution of the helicity and the total angular momentum carried by quarks, starting at μ = 2 GeV with input values from lattice calculations and evolving to low scales in order to compare with quark models. We note that significant deviations between LO and higher orders (NLO and NNLO) were found below μ ∼ 500 MeV, and we anticipate that we will find the same in the unpolarised sector studied here.
The present work is structured as follows. In Sect. 2, we discuss some obvious caveats of our study. We then introduce our notation and briefly recall the evolution of Mellin moments in Sect. 3. In Sect. 4, we present the different PDF sets that we use to compute the starting values for evolution. In the same section, we quantify to which parton momentum fractions x one is most sensitive in a given Mellin moment, and we take a brief look at the running of the strong coupling. In Sect. 5, we investigate in detail the evolution of Mellin moments to low scales. The main findings of our study are summarised in Sect. 6.

Caveats
Our study is subject to several caveats, which we now briefly discuss. Perhaps the most obvious one is our use of perturbation theory down to rather low scale. In fact, there is a substantial body of work on the behaviour of α s in the low-scale regime; see, for instance, the review [30] and the recent lattice studies [31][32][33]. However, to the best of our knowledge, there are no corresponding studies for the scale dependence of PDFs or their Mellin moments. To assume that PDFs evolve as given by perturbation theory in a region where the running of α s (μ) is strongly affected by non-perturbative effects would in our view require some motivation. In the present work, we use perturbation theory to evolve both the strong coupling and Mellin moments to low scales with the aim to see what one obtains in such a scenario, without strong claims that this is a valid approximation at a given scale μ. In doing so, we follow the procedure adopted in many works that connect quark models at low scale with PDFs (see the papers cited in the introduction). By comparing perturbative results at three different orders, we will in fact get some indication of how stable the perturbative expansion is at a given scale for the quantities we are interested in.
Another caveat concerns the renormalisation scheme used to define PDFs and α s . We will use the MS scheme throughout this work. This choice is dictated by practical considerations, as it is in this scheme that the DGLAP splitting functions are available up to NNLO. The strength of this scheme is its suitability for higher-order perturbative computations, but as a downside it does not readily offer an intuitive interpretation of renormalised quantities. One might think of other schemes for defining PDFs, e.g. the DIS scheme, but this will not be pursued in the present work. (Notice that, whilst quark and antiquark densities have a rather straightforward physical meaning in that scheme, the same is not true for the gluon distribution.) We note that PDFs and the running coupling at LO play a special role in this context, since they are the same in a large class of renormalisation schemes. Returning to the discussion in the previous paragraph, one might ideally want to define a scheme that allows for a physical interpretation of PDFs and that can be used in a non-perturbative setting, but such an endeavour is well beyond the scope of this work.
Finally, the very idea that antiquark or gluon distributions in a hadron are zero at a certain scale μ c is somewhat problematic regarding their physical interpretation. This is because these distributions will then typically become negative at scales just below μ c , at least in some x range. We will see this happen at the level of their Mellin moments. The interpretation of PDFs as probability densities is then lost at scales below μ c . One may mitigate this problem by considering a scale just above μ c , where antiquark or gluon densities would be nonzero but small (compared with their values at higher scales, or compared with the quark distributions at the same scale).

Evolution of Mellin moments
Let us briefly recall how the Mellin moments of PDFs depend on the renormalisation scale. Throughout this work, we consider distributions for n f = 3 active quark flavours and limit our attention to the flavour singlet combinations of quark and antiquark distributions. One could extend our study to individual quark flavours, investigating, for instance, whether moments of the distribution s(x, μ) +s(x, μ) vanish at some low scale. In such a case, strangeness in the proton would be generated by perturbative evolution. To pursue this is, however, beyond the scope of this paper.
We define Mellin moments for flavour summed quark and antiquark distributions, and for the gluon distribution, recalling that all quantities are renormalised in the MS scheme. The valence combination Q − Q evolves as where the non-singlet anomalous dimension γ ns ( j, μ) depends on μ via the scale of α s . In the singlet sector, we have the matrix equation The anomalous dimensions have a perturbative expansion, which reads Using the renormalisation group equation dα s /d log μ 2 = β(α s ), we rewrite (3) and (4) as evolution equations in the variable α s , which gives and an analogous equation for the singlet sector. We use the perturbative expansion of β(α s ) at the same order as for the anomalous dimensions and solve the equations numerically using the classical four-step Runge-Kutta method. For brevity, we suppress the dependence of the Mellin moments on μ (or on the corresponding value of α s ) henceforth. The non-singlet equation (6) can of course be solved analytically by straightforward integration. The same holds for the singlet equation if j = 2, in which case G(2) can be eliminated by using the momentum sum rule Q(2)+ Q(2)+ G(2) = 1. We use these analytic solutions as cross-checks of the numerical ones.
At odd values of j, the combination Q( j) − Q( j) of moments is connected with the proton matrix elements of local twist-two operators, and the same holds for Q( j) + Q( j) and for G( j) if j is even. In our study, we are, however, interested in the antiquark moments Q( j) by themselves, which are not connected with local twist-two operators at any value of j. Furthermore, for reasons discussed in Sect. 4.1, we will also consider non-integer values of j. We obtain the anomalous dimensions needed in (3) and (4) by computing the relevant Mellin moments of the DGLAP splitting functions. For the splitting functions up to NLO (k = 1), we use the exact expressions given in [34], whereas for the NNLO (k = 2) splitting functions we take the parametrisations given in [35,36], which approximate the exact kernels and have a much simpler functional form. In both cases, the Mellin moments are easy to compute. We cross-checked our results for the j = 2 anomalous dimensions in the singlet sector against the analytic expressions given in [37]. For the NNLO coefficients, we also verified the constraints from fermion number and momentum conservation for our numerical results and find them to be satisfied within better than 2 × 10 −3 for n f = 3. We take this as evidence that the approximate forms of the NNLO splitting functions in [35,36] are sufficiently accurate for our purposes. In Table 1, we give the numerical values of the perturbative expansion coefficients in (5) as functions of n f . We note that approximate expressions for the DGLAP splitting functions in the non-singlet sector are available at N 3 LO [38]. Furthermore, the N 3 LO coefficients γ (3) i ( j) in the singlet sector have been given in [39] for j = 2 and n f = 4. This is, however, not sufficient for extending the expansion coefficients given in Table 1 to the next order, k = 3, and we limit our present analysis to k ≤ 2.

Parton densities, their moments and the running coupling
We perform our study with a wide range of current PDF sets, which should reflect the current knowledge and uncertainties of unpolarised parton densities. 1 In Table 2, we list these sets together with some of their characteristic parameters. In the first column, we give the full name of a set in the LHAPDF library [40], from which we take the numerical values of all PDFs. The second column shows the "short names" that will be used to refer to a given set throughout this paper. A number of comments are in order.
• The 2018 Review of Particle Physics [1] gives α s (M Z ) = 0.1181 ± 0.0011 as world average for the strong coupling at the Z mass, which corresponds to values from 0.1159 to 0.1203 at 2σ accuracy. To assess the impact of the α s value within a single PDF fitting approach, we take the NNLO sets of NNPDF [41] for α s (M Z ) = 0.116, 0.118, 0.120. We note that the value α s (M Z ) = 0.11471 in the ABMP set [42] is even smaller. • As starting scale for evolution, we take μ 0 = 1.3 GeV for all PDF sets. At this scale, all sets using a variable flavour number scheme have n f = 3 active quarks. The charm quark mass m c in the ABMP set is smaller than μ 0 , but this set uses the fixed flavour number scheme with n f = 3.  • The default PDF set in the NNPDF study [41] has an intrinsic charm distribution and is hence not available for n f = 3. Evolving a n f = 4 set down to low scales makes little physical sense. We therefore take the "perturbative charm" variant of that study, in which charm distributions are generated by evolution, as is the case for all other PDF sets in our study.

Table 2
Overview of the PDF sets considered in the present study. It is understood that the number of active quark flavours is n x min ABMP16_3_nnlo ABMP [42] 0.340    Fig. 1 Fractional contribution of different x intervals to selected Mellin moments. In some cases, the very small x or the large x region gives a negative contribution, indicating that the corresponding PDFs cannot be interpreted as number densities in that region. The PDF sets are grouped according to their perturbative order (LO 1-4, NLO 5-10,  • In the tradition of PDF fits by the Dortmund group, the JR study [8] takes an initial condition at very low scale, namely at Q 2 0 = 0.8 GeV 2 . Such an approach is of obvious interest in the context of our study. To assess its impact on the moments of PDFs, we also include the set with a more conventional starting scale of Q 2 0 = 2.0 GeV 2 from the same study. We refer to the respective sets as "JR 08" and "JR 20".
• We include PDF sets from CJ [43], because that study pays particular attention to the region of large x. We find that the PDF error bands for these PDFs, as given by the LHAPDF interface, are considerably smaller than those of any other set we studied. We have not investigated the reasons for this and decided not to show the corresponding sets in our plots. We will, however, include the CJ15 sets in our discussion later on. • The LO set of CT [44], as implemented in LHAPDF, does not give any PDF uncertainties, and we therefore exclude it from our study.

PDF moments and x ranges
We consider a collection of moments, from j = 1.5 to j = 3 in steps of 0.5, in order to be sensitive to PDFs in a wide range of momentum fractions x. To quantify this sensitivity, we divide the x range into four intervals, to which we, respectively, refer as "very small x region", "small x region", "valence region", and "large x region" in the sequel. We then take the PDFs at the starting scale μ 0 of our study and determine the contribution of these different x intervals to each Mellin moment. The result is shown for selected moments in Fig. 1 and summarised in Table 3. Notice that the gluon moments have the strongest variation between different sets, with less variation for antiquarks and even less for quarks. We see that G(2) and Q(2) receive comparable contributions from the small x region and the valence region. To have quantities that are dominated by the small x region, we resort to non-integer moments with j = 1.5. (The next smallest integer j = 1 is not an option, because the corresponding Mellin integrals are in general divergent for G,Q, and Q.) We note that the very small x region, in which fitted PDFs are essentially unconstrained by data, plays only a minor role in all moments we consider. The high moments with j = 2.5 and j = 3 are increasingly dominated by the valence region. The region of large x, where gluon and antiquark distributions are again poorly known, plays only a minor role in the G and Q moments. To summarise, we find that at our starting scale μ 0 = 1.3 GeV, the Mellin moments we consider offer a reasonably differential sensitivity to the region 10 −4 < x < 0.5, in which PDFs are reasonably well known. Moments with lower or higher j will be more strongly affected by PDF uncertainties. When going to lower scales, one should keep in mind that the PDFs are in general shifted towards higher x values by backward evolution. The relative importance of the different x regions for a given j will then change.
We note that for some PDF sets, the gluon distribution at μ 0 = 1.3 GeV becomes negative at very small x. Some sets have a zero crossing of g(x, μ 0 ) at x below 10 −4 . For other sets (such as HERAPDF NNLO and MMHT NNLO), the zero crossing occurs already for x below 10 −3 , and we see in Fig. 1 that the very small x region gives a negative contribution to the lowest moment G(1.5). Such a behaviour is a clear indicator that the corresponding PDFs no longer admit a probability interpretation in the corresponding region.
When computing Mellin moments, we truncate the integral over x at the smallest value x min for which LHAPDF provides PDF values for a given set, thus avoiding an extrapolation of parton densities down to x = 0. The values of x min are given in the last column of Table 2 and range from 10 −9 to 10 −6 . To obtain a conservative estimate of the uncertainty on the moments due to this truncation, we recompute them with a lower integration limit of 10 x min instead of x min . The resulting change is less than 0.1% for all moments with j ≥ 2. For the j = 1.5 moments, it is less than 3% with the following exceptions. The antiquark moment Q(1.5) of the NNPDF LO set changes by 4%, which is negligible compared with the huge error on this moment due to the PDF uncertainty at x > x min . The gluon moment G(1.5) changes by 5% for HERAPDF NNLO and by 4%, 8%, and 7% for the MMHT sets at LO, NLO, and NNLO, respectively. Our conclusions in Sect. 5 are not affected by these somewhat larger truncation uncertainties.
When evolving the Mellin moments to lower scales, we use anomalous dimensions at fixed order in perturbation theory. One may wonder whether some type of all-order resummation would improve perturbative convergence. We argue that this is not the case: small-x logarithms α s log(x) correspond to powers of α s /j in Mellin space and hence do not require resummation for j > 1. Large-x logarithms, which correspond to powers of α s log 2 j [47,48], do not appear in anomalous dimensions for the evolution of parton densities [35,36,49], and in any event, log j is not large for j ≤ 3.

The running coupling
Let us take a brief look at the running coupling α s (μ). In Fig. 2a, we show the evolution of α s (μ) down to low scales at different orders in the perturbative expansion of β(α s ), which is available up to five-loop order [50], i.e. up to N 4 LO. For definiteness, we take in this plot a common value α s (1.3 GeV) = 0.378 for all orders. (This is the value of the CT NNLO set as seen in Table 2.) The plot looks very similar if instead one takes a common value α s (M Z ) = 0.118 for all orders and sequentially evolves and matches the coupling from n f = 5 down to n f = 3. We observe that down to about μ ∼ 0.7 GeV, the difference between different orders is small and decreases with the order, indicating satisfactory convergence of the perturbative expansion. For lower scales, however, the different orders differ more and more strongly.
In Fig. 2b, we show the running of α s (μ) at NNLO with starting values at μ 0 = 1.3 GeV corresponding to those in different PDF sets. We see that the moderate spread in α s values at μ 0 rapidly increases when evolving to lower scales.

Evolution to low scales
We have now everything in place to investigate the evolution of Mellin moments for different j from their starting values at μ 0 = 1.3 GeV down to small scales. We focus on the moments  Table 2, the curves for the remaining NNLO sets of our study are in between those shown in the plot G( j) and Q( j). As a function of μ, PDF moments exhibit a very steep behaviour at scales where α s (μ) starts to diverge, as is seen in Fig. 6. To display the low-scale behaviour of the moments in a clearer way, we will in general plot them as functions of α s instead of μ.
We remark that the lowest moments G(1.5) and Q(1.5) are consistent with zero within huge errors for three sets, namely for NNPDF LO, CT NLO, and CT NNLO. This already holds at μ 0 and reflects the huge uncertainties of the antiquark and gluon distributions at low x in these sets. Evolving G(1.5) and Q(1.5) to yet lower scales then gives no further useful information. For simplicity, we will not explicitly mention these cases in the following discussion.

Comparison of different orders
If we study PDF moments as functions of α s , then the quantity that drives their evolution is γ i (α s )/β(α s ) according to Eq. (6) and its analogue for the singlet channel. In Fig. 3, we show this quantity for the three perturbative orders considered in our study, multiplied with an additional power of α s so that the LO curves are constants. The channels and moment indices j shown in the figure are representative of the wide range of patterns we observe. In some cases, there are only moderate differences between different orders, in other cases the difference between LO and NNLO amounts to a factor around two at the highest α s shown in the figure, and in yet other cases one finds that the NNLO curve changes sign at some intermediate value of α s . Let us emphasise that by showing the curves up to α s = 3, we do not mean to imply that perturbation theory is valid up to that value. Indeed, both Figs. 2 and 3 suggest that for the quantities we are studying here, perturbation theory breaks down well before that value of α s is reached.
To assess the difference in the evolution behaviour of Mellin moments at different orders, we focus on those sets that provide PDFs at all three orders, i.e. on HERAPDF, MMHT and NNPDF (see Table 2). A selection of moments for the HERAPDF sets is shown in Fig. 4. For G( j) with j ≥ 2, we find the NLO and NNLO moments to be rather close to each other, whilst the LO moments are clearly different from these at large α s . For G (1.5), all orders differ quite noticeably. The situation is similar for MMHT and NNPDF. For the antiquark moments, we find in general larger differences between different orders, with a pattern that depends on j and also on the considered PDF set. This stronger dependence on the initial conditions is perhaps not surprising, given that Q is the difference between the combinations (Q + Q)/2 and (Q − Q)/2, which evolve independently of each other. We note that the width of the error bands for the moments increases with α s , but in most cases it does so rather slowly. The same is observed for the other PDF sets and is consistent with the numerical stability of backward evolution for Mellin moments.

Comparison of different PDF sets
We now discuss how given Mellin moments compare between different sets. In addition to the parametric PDF errors provided by each set, the difference between sets may be taken as an indication of how well a given Mellin moment is known. The situation is shown in Fig. 5 for j = 2 at LO and in Fig. 6 for j = 3 at all orders. In this case, we plot moments against μ rather than α s , because different PDFs sets have different values of α s (μ) at a given scale μ. Starting with the LO sets, we find that G( j) evolves to a zero value and then becomes negative for j = 1.5 and 2. For j = 2.5, it either increases or decreases with α s , depending on the set, and for j = 3 it increases or remains flat. The moments Q( j) with j = 1.5 and 2 increase with α s . For j = 2.5 and 3, the antiquark moments evolve to zero for MMHT, whereas for HERAPDF and NNPDF they decrease only slightly or remain flat. The behaviour of Q (2) and G(3) in Figs. 5b and 6a alone implies that there is no LO set in which either the gluon or the flavour sum of antiquark distributions evolves to zero at some low scale. This remains true if we include the CJ LO set in our considerations. (See our remark on that set in Sect. 4.) Moving on to the NLO sets, we find that G( j) has a zero crossing for j = 1.5 and 2 in all sets, whereas for j = 2.5 and j = 3 this only happens for some of the sets. As was the case for LO, the moments Q( j) with j = 1.5 and 2 increase with α s . We find no zero crossings for Q (2.5), whereas in some PDF sets the error bands of Q(3) just touch zero at the lowest scales, as seen in Fig. 6d. None of the NLO sets is hence compatible with the flavour sum of antiquark distributions vanishing at low scales.
Turning to the NNLO sets, we observe that in all sets all moments G( j) decrease with α s , either crossing zero, or coming close to zero, or being consistent with zero within their error bands. The antiquark moments with j = 1.5, 2, and 2.5 decrease with α s , with some of them going down to zero and others not. For the moment Q(3), only HERAPDF is consistent with zero at very low scale, whilst the other sets are not, as seen in Fig. 6f. Whether any NNLO set admits the gluon or antiquark distributions to evolve to zero is investigated more closely in the next subsection.

Vanishing gluon or antiquark moments
If g(x, μ) or the sumū(x, μ) +d(x, μ) +s(x, μ) is zero at some scale μ, then all Mellin moments of these distributions must vanish at that scale, or equivalently at the corresponding value of α s (μ). In order to test this hypothesis, we compare the different moments G( j, α s ) and Q( j, α s ) for each PDF set in turn.
Starting with the gluon, we see in Fig. 7 that for ABMP, JR NNLO 08, and HERAPDF NNLO, all moments G( j) have a zero crossing but are never consistent with zero at the same value of α s . In particular, the lowest moment G(1.5) has already evolved to negative values at α s values where the higher moments become zero, so that g(x) must be negative in some  Fig. 7d for the 118 set. In all other PDF sets (including CJ), we find that G(3) stays positive. (For MMHT NNLO, the error band just touches zero at α s ≈ 3.) We thus find no set that is compatible with a vanishing (or very small) gluon distribution at low

Conclusions
In this work, we investigate the possibility that the gluon or the antiquark PDFs in the proton go to zero when evolved to low scales, using the DGLAP equations with splitting functions computed in perturbation theory. To do so, we take a variety of current PDF sets and compute Mellin moments of PDFs for j = 1.5, 2, 2.5, and 3 at a reference scale μ 0 = 1.3 GeV. We then evolve these moments down to lower scales, stopping at the point where α s (μ) = 3. Our comparison of these moments at LO, NLO, and NNLO indicates that perturbation theory ceases to converge at much lower values of α s , so that considering even larger values would make little sense.
In several NLO and NNLO PDF sets, the gluon moments for all considered j go to zero under evolution, but they do so at different scales μ. No PDF set is found to be compatible with the gluon PDF vanishing at any scale. We note that the j = 1.5 moment eventually evolves to negative values in almost all sets at some scale μ below 0.8 GeV, which indicates that g(x) < 0 at low x and hence cannot be interpreted as a probability density.
For antiquarks, we take the flavour sumū(x) +d(x) +s(x) and find that at least some of its Mellin moments remain positive up to the largest α s we consider. An exception to this is the HERAPDF NNLO set, for which all antiquark moments congregate around zero for α s 2.8. At these scales, however, the gluon density in the same set has negative Mellin moments, so that we cannot interpret the PDFs as densities.
Our findings are fully consistent with those of the Dortmund group [5][6][7][8], who concluded from their fits that PDFs that describe high-energy scattering data cannot be generated by perturbative radiation from an input that involves only valence quark densities at some low momentum scale. We think that our study is a valuable complement to the Dortmund approach, in particular, because we use the results of a broad array of PDF fits, which differ in their data selection, parametrisation of PDFs, as well as the details of the theory description of hard cross sections.
Which scenarios does this situation leave for connecting the PDFs determined from data at high scales with PDFs computed in models at low scales? We see several obvious possibilities. The least dramatic one would be to consider PDFs defined in a perturbative renormalisation scheme other than MS and to identify these PDFs with the results of low-energy computations. Perhaps more plausible is that the evolution of PDFs is modified by non-perturbative effects at low scales, as has been argued for the case of the running coupling [30]. This is also suggested by our comparison of Mellin moments evolved at different orders, and by the study [28] of evolution in the polarised parton sector. Finally, it may well be that even at the lowest scales one can sensibly consider, the parton content of the proton is not limited to just "valence quarks" but involves antiquarks or gluons or both, as several low-energy models suggest.

Note added in proof
After this work was completed, we were made aware of the paper [51]. Using PDF sets available at the time, that study evolved PDFs in x space down to low scales. In particular, parton distributions were found to become negative below μ ∼ 600 MeV. The qualitative conclusions of [51] agree with those of our work, disfavouring the idea that PDFs at high scales can be obtained from DGLAP evolution of PDFs computed in quark models.