Leading order DGLAP evolution Consider the leading-order (LO) DGLAP ((Dokshitzer–Gribov–Lipatov–Altarelli–Parisi) equation
$$x \mu^2 \frac{d xg(x,\mu^2)}{d\mu^2}= \alpha_s \int_x^1 dz P_{gg}(z) \frac{x}{z} g(\frac{x}{z}, \mu^2) + \dots\tag1$$ 
Define an unintegrated parton density function (PDF) by
$$F(x,\mu^2) \equiv x \frac{dg(x,\mu^2)}{d\mu^2},$$
with the initial condition
$$\mu^2 F^{(0)}(x,\mu^2) = \theta(1-x) \theta \left(\frac{\mu^2}{Q_o^2}-1\right),\tag2$$
where $Q_o$ is some non perturbative scale. The first order contribution to $F$ is then given by inserting $(2)$ into $(1)$ so that
$$\mu^2 F^{(1)}(x,\mu^2) \approx \int_x^1 dz P_{gg}(z) \int_{Q_o^2}^{\mu^2} dk^2 \alpha_s(k^2) F^{(0)}(\frac{x}{z}, k^2)$$
My questions are 
1)From the definition of $F$ it follows that $$g(x,\mu^2) = x^{-1} \int^{\mu^2} dk^2 F(x,k^2)$$ so I can see how the integral in the last display arises but what I don't see is why does this correspond to the first order contribution to $F$? Doesn't equation (1) tell us that the evolution of the LO gluon PDF is given by the integral convolution of the LO splitting function kernels with the LO gluon PDF itself? So where does next-leading order (NLO) effects, i.e that suggested by 'first order in $F$' arise?
2)What are the physical meanings of the $F^{(i)}$? Are these defined through some expansion and, if so, what is this expansion in?
 A: First, the DLAP equation at any order can be written
$$\mu^2x\frac{dg(x,\mu^2)}{d\mu^2} = \int_x^1 dz P_{gg}(z) \frac{x}{z} g\left(\frac{x}{z}, \mu^2\right) + \cdots$$
Then, the kernel $P_{gg}(z)$ is expanded in powers of $\alpha_S$:
$$P_{gg}(z) = P_{gg}^{(0)}(z)\frac{\alpha_S}{2\pi} + P_{gg}^{(1)}(z)\left(\frac{\alpha_S}{2\pi}\right)^2 + \cdots$$
Then LO means only the first term and NLO means the first two terms are included. Thus I think that the wording "first order contribution" used in the work you study for $F^{(1)}$ is a different notion. Basically, this is a scheme to solve the DGLAP equations by successive approximations: the idea would be that $F^{(n)}$ converges to the solution as $n$ increases, I think.
This should work as follow. After defining (beware I added a factor $\mu^2$ on the right-hand side so as to keep the equations more symmetric later on)
$$F(x,\mu^2) \equiv \mu^2 x \frac{dg(x,\mu^2)}{d\mu^2},\label{defF}\tag{II}$$
the DGLAP equation reads
$$F(x,\mu^2)=\int_x^1 dz P_{gg}(z) \frac{x}{z} g\left(\frac{x}{z}, \mu^2\right) + \cdots\label{dglapalt}\tag{I}$$
we start from a first approximation given by your equation (2), to which correspond an approximation $g^{(0)}$ of $g$ given by
$$xg^{(0)}(x,\mu^2) = \int_{Q_0^2}^{\mu^2} \frac{dk_0^2}{k_0^2} F^{(0)}(x, k_0^2).$$
I simply integrated with respect to $\mu^2$. This approximation can then be inserted in the right-hand side of the DGLAP equation (\ref{dglapalt}) to give 
$$F^{(1)}(x,\mu^2) \approx \int_x^1 dz_0 P_{gg}(z_0) \int_{Q_o^2}^{\mu^2} \frac{dk_0^2}{k_0^2} F^{(0)}\left(\frac{x}{z_0}, k_0^2\right)+\cdots$$
Then from $F^{(1)}$, we can get $g^{(1)}$ by integrating with respect to $\mu^2$ again, 
$$xg^{(1)}(x,\mu^2) = \int_{Q_0^2}^{\mu^2} \frac{dk_1^2}{k_1^2} F^{(1)}(x, k_1^2)=\int_x^1 dz_0 P_{gg}(z_0)\int_{Q_0^2}^{\mu^2}\frac{dk_1^2}{k_1^2}\int_{Q_0^2}^{k_1^2} \frac{dk_0^2}{k_0^2} F^{(0)}\left(\frac{x}{z_0},k_0^2\right) + \cdots$$
then insert $g^{(1)}$ in the right-hand side of (\ref{dglapalt}) to get you $F^{(2)}$,
$$\begin{aligned}
F^{(2)}(x,\mu^2)&=\int_x^1 dz_1 P_{gg}(z_1) \frac{x}{z_1} g^{(1)}\left(\frac{x}{z_1}, \mu^2\right) + \cdots\\
&=\int_x^1 dz_1 \int_{z_1}^1 dz_0 P_{gg}(z_1) P_{gg}(z_0)\int_{Q_0^2}^{\mu^2}\frac{dk_1^2}{k_1^2}\int_{Q_0^2}^{k_1^2} \frac{dk_0^2}{k_0^2} F^{(0)}\left(\frac{z_1}{z_0},k_0^2\right) + \cdots
\end{aligned}$$
We can then continue by computing $g^{(2)}$, etc. Each step of this procedure add one more splitting. This can probably be connected to the origin of the DGLAP equations, which appear as the result of the resummation of a class of amplitudes featuring 1, 2, 3, etc colinear gluon emissions.
