Are parton distribution functions non-negative? Let $f_i(x)$ be a parton distribution function as known from QCD factorisation theorems. Is $f_i(x)$ non-negative for $0<x<1$? If so, how can this be seen from the definition of PDFs in terms of operator products?
 A: By definition parton distribution function is a re-scaled probability density :
$$ f_i(x) = K\, |\phi |^{2} $$
So YES, it must be non-negative, however it can be greater than $1$.
Also take a look into this research :

where authors claim that :

The VFN scheme is valid only at asymptotically large values of Q and
cannot be routinely extrapolated to the low-Q region.

So it can't take negatives values for sure.
A: OK, as the question eventually morphed, the answer is


*

*For plain (non phase space GPDs!) pdfs, they really are scaled probability distribution functions, always positive semidefinite, and their integrals over x are positive numbers.


They are an assumed input (not an output) of parton estimates, fitted to data. 
Typically, as, e.g., in Schwartz 32.1.4, they count gluons and constituent and sea quarks. You know there is a constituent d in the proton. If there were no sea, you'd have 
$$
\int dx ~~f_d(x) =1,
$$
and the function would be a bona-fide probability distribution: it would map what fraction of the proton's momentum that quark is likely to carry. 
But there is a sea, so also an extra contribution to $f_d(x)$, which actually blows up at small xs. However, that extra contribution is exactly matched by $f_{\bar d} (x)$, so that $\int dx ~  f_d(x) =1+ \int dx  ~f_{\bar d}(x)$,
$$
\int dx ~ (  f_d(x)  -f_{\bar d}(x))=1,
$$
implying that the actual normalization of the probability distribution $f_d(x)$ is a positive number, left undisclosed, and dependent on other features of the problem (Q). (If you had it, you could normalize your $f_d(x)$ with it, assuming that mattered to you.) 
Likewise, now, since there are 2 valence u quarks, it makes sense to scale their probability distribution normalization by a factor of 2, so you count them in all, $\int dx ~f_u(x)=2$, but, given the sea, again,
$$
\int dx ~ (  f_u(x)  -f_{\bar u}(x))=2,
$$
etc; so that, at the end of the day,
$$
\int dx ~ (   f_u(x)  -f_{\bar u}(x)  +    f_d(x)  -f_{\bar d}(x) +f_s(x)  -f_{\bar s}(x)+  f_c(x)  -f_{\bar c}(x)+ ... )                    =3,
$$
for 3 constituent quarks in all. 


*

*Each f is a positive probability, suitably normalized, by assumption, and obeys the rules of probability.


[By contrast, GPDs in phase space are not, and violate probability axioms with aplomb, but since you did not ask, let's leave them out for now...]
A: In light cone gauge parton distribution functions are matrix elements of number operators inside a hadron, so they appear to be manifestly positive, as expected from the naive parton model. There is a subtlety, because these matrix elements are divergent, and require counter terms. As a result, I do not think that there is a formal proof of positivity. 
Of course, parton distributions determine cross sections, and cross sections must be positive. 
