How to derive operator form of the parton distribution function A similar question is found here in Stackexchange a year ago without any response.
I am following the formulation of the parton densities from the handbag diagram in Collins Handbook of Perturbative QCD, Chapter 6. I know that the quark density can be written as (left as an exercise 6.2)
$f_{j/h}(\xi) = \int\frac{dw^-}{2\pi}e^{-i\xi P^+w^-}
<P|\bar{\psi}_j(0, w^-,\textbf{0}_T)\frac{\gamma^+}{2}\psi_j(0)|P>$ (6.31).
I am having difficulty figuring out where $\gamma^+$ comes from. 
Also in Chapter 7, 7.5.4 gauge-invariant quark density in QCD, they make make the (bare) quark density function gauge invariant by inserting a Wilson line in between the two fermion fields:
$f_{(0)f/h}(\xi)=\int\frac{dw^-}{2\pi}e^{-i\xi P^+w^-}
<P|\bar{\psi}_j(0, w^-,\textbf{0}_T)\frac{\gamma^+}{2}W(w^-,0)\psi_j(0)|P>$ (7.40).
Where exactly does the Wilson Line come from? I feel uncomfortable just saying that the Wilson Line is put there by hand to make that matrix element gauge-invariant, and would like to see a derivation. 
 A: There is a pedagogical derivation of the operator form of the parton distribution function presented in e.g Schwartz involving a gauge invariant matrix element with light cone separated quark fields but again with the ad hoc insertion of a Wilsonian line connecting them.
From a physically motivated perspective within the framework of collinear factorisation, Schwartz finds $$f(\xi) = \int \frac{\text{d}\lambda}{2 \pi} e^{- \mathrm{i} \lambda \xi n \cdot P} \langle P | \bar{\psi}(\lambda n) \gamma^0 \psi(0) | P\rangle $$
One can write a generic $4$ dimensional vector in a Sudakov basis spanned by vectors e.g $n \sim (\Lambda,0,0,\Lambda)$ and $\bar n \sim (\Lambda,0,0,-\Lambda)$ such that $$A^{\mu} = (n \cdot A) \bar n^{\mu} + (\bar n \cdot A) n^{\mu} + A^{\mu}_{\perp}.$$  Setting $\Lambda$ as unity we make the identification of the Sudakov basis with the light cone basis as used in Collins’ formulation.  That is to say, $A^+ \equiv (n \cdot A)$ and $A^- \equiv (\bar n \cdot A)$ so that $$A^{\mu} = A^+ \bar n^{\mu} + A^- n^{\mu} + A^{\mu}_{\perp}.$$
In an analogous fashion, one can do the same with a gamma matrix (i.e send $A \rightarrow \gamma$)  and with an appropriate normalisation of the Sudakov vectors may identify $$\gamma^0  = \frac{1}{\sqrt{2}} \left(\gamma^+ + \gamma^-\right).$$ This of course also follows from the canonical light cone coordinate definitions: $$\begin{cases} \gamma^+ = \frac{1}{\sqrt{2}} \left(\gamma^0 + \gamma^3 \right) \\ \gamma^- = \frac{1}{\sqrt{2}} \left(\gamma^0 - \gamma^3\right) \end{cases}$$
As per Schwartz, if the proton is going mostly in the $\bar n$ direction, then the projection of $\psi$ in the $n$ direction is suppressed (contributing at higher twist outwith the famous hard and soft sector factorisation). This all means to say $(\bar n \cdot \gamma) \psi = \gamma^- \psi \approx 0$  and we obtain $$\gamma^0 \psi \approx \gamma^+ \psi$$
A: As I am currently reading Collins book and stumbled upon this post I give my understanding of how the $\gamma^+$ arises.
Collins explaines in previous sections that the leading regions of the whole Feynman diagram comes from a subspace of loop momentum space, where the lines in the lower bubble of handbag diagram are approximately collinear to the target hadron, i.e. they carry momenta of order $P \sim (Q,m^2/Q,m)$,
where $m$ is a generic mass scale of the particles which is small compared to $Q$.
Let $L$ be the amplitude for the lower subgraph, which we can decompose as:
$$
L = A + \gamma_5 B + \gamma_{\mu} C^{\mu} + \gamma_{\mu} \gamma_5 D^{\mu} + \sigma_{\mu \nu} E^{\mu \nu}
$$
Collins' all order argument is that of 6.1.2. In the rest frame of the subgraph, i.e. when the hadron momentum scales like $P \sim (m,m,m)$, every coefficient is of the same size. When we now boost to the our original frame (the Breit frame), where the hadron has collinear momentum, we greatly increase the $+$ component and decrease the $-$ component. A vector  $v^{\mu}$ transforms under this boost as
$$
(v^+,v^-,\mathbf{v}_T) \rightarrow (v^+ e^y,v^- e^{-y},\mathbf{v}_T) \sim (v^+ Q/m,v^- m/Q,\mathbf{v}_T)
$$
where $y = \frac{1}{2} \ln \frac{P^+}{P^-} \sim \ln Q/m $ is the rapidity of $P$. The statement is now that the coefficients $A,B,C^{\mu}$,etc. transform according to their Lorentz group representation, $A$ as scalar, $C^{\mu}$ as vector, etc. So the leading contributions will come from the terms $C^+, D^+, E^{+j}$ ($j=1,2$). In particular the term proportional to $\gamma^+$ can be obtained from $L$ by
$$
C^+ = \frac{1}{4} \text{Tr } \gamma^+ L
$$
The other contributions $D^+, E^{+j}$ correspond to the polarized PDFs.
This is actually almost as detailled as Collins gets in his book, at least up to the point where I have read so far. A similar boost argument is also used in the derivation of the power counting laws which are very important for Factorization proofs in pQCD. If still find this quite handwavy (similar to the argument in Schwartz's book) so I would be interested if someone knows any source where this is shown more rigorously.
Also note that the example in chapter 6 is for a super-renormalizable non-gauge theory and it does not necessarily work exactly the same in QCD.
As for the Wilson line: I know that it "sums up" the interactions of collinear gluons connecting the upper hard subgraph and the lower collinear subgraph, so the justification of its appearance goes well beyond just gauge invariance! The diagrams with these collinear gluon lines are not power-suppressed (as shown in section 5.8) and should be regarded in the leading order factorization formula. It is discussed in the latter sections of the book and it is probably to complicated to expand on this here.
