How are DGLAP splitting kernels probabilites for parton branchings? In the literature on parton evolution one often reads that the DGLAP splitting kernels $P_{ij}$ represent the probabilities for a parton of type $j$ to split into a parton of type $i$ and another parton.
At leading order the kernels are:
\begin{align}
P_{qq}(z)&=C_F\left(\frac{1+z^2}{[1-z]_+}+\frac{3}{2}\delta(1-z)\right),\\
P_{qg}(z)&=T_F\left(z^2+(1-z)^2\right),\\
P_{gq}(z)&=C_F\frac{1+(1-z)^2}{z},\\
P_{gg}(z)&=2C_A\left(\frac{z}{[1-z]_+}+\frac{1-z}{z}+z(1-z)\right)+\frac{\beta_0}{2}\delta(1-z),
\end{align}
where $\beta_0=\frac{11}{3}C_A-\frac{4}{3}T_FN$ (see pp.681) in M.D. Schwartz: Quantum Field Theory and the Standard Model.
If the splitting kernels are to be interpreted as the probability densities for a parton of type $i$ to branch into a parton of type $j$, we should have $\int^1_0 dz \sum_j P_{ji}(z)=1-P(\mathrm{no branching})$.
But integrating the above expressions gives
\begin{align}
\int^1_0dzP_{qq}(z)&=0,\\
\int^1_0dzP_{qg}(z)&=\frac{2}{3}T_F,\\
\int^1_0dzP_{gq}(z)&=C_F\left(2\int^1_0\frac{dz}{z}-\frac{3}{2}\right),\\
\int^1_0dzP_{gg}(z)&=C_A\left(2\int^1_0\frac{dz}{z}-\frac{11}{6}\right)-\frac{2}{3}NT_F,
\end{align}
i.e. two of the kernels are clearly not integrable over $[0,1]$, which is a prerequisite for any function to be a probability density.
How can the kernels still be interpreted as probability densities? Or are they maybe no probability densities to begin with? Maybe they can only be interpreted as such, if they are defined on a domain $[x,1]$ with $0<x<1$?
 A: The splitting functions tell you the probability for a given radiance to occur, and some exhibit manifest divergences (e.g. as $z \rightarrow 0$) signifying a growth in that particular branching in a certain part of the phase space. The correct statements are $$\int_0^1 z \, \text{d}z \, (P_{q \leftarrow q}(z) + P_{g \leftarrow q}(z)) = 0\,\,\,\,\,\,\,;\,\,\,\,\,\,\,\,\int_0^1 z\, \text{d}z\, (2n_f P_{q \leftarrow g}(z) + P_{g \leftarrow g}(z)) = 0 $$
They may be derived by totalling the sum of all  fractional constituent nucleon momentum to unity, means $$ \int_0^1 \text{d}x \,x \, \left[ \sum_{i} \left( q_i(x) + \bar{q}_i(x) \right) + g(x) \right] = 1 $$ and taking the derivative w.r.t $\mu^2 \,d/d \mu^2$  in each of the terms as follows $$\sum_i \mu^2 \frac{d}{d \mu^2} ( q_i(x) + \bar{q}_i(x) ) = \frac{\alpha_s}{2 \pi} \sum_i \left( q_i \otimes P_{q \leftarrow q} + 2 g \otimes P_{q \leftarrow g} \right),$$ $$ \mu^2 \frac{d}{d \mu^2} g(x) = \frac{\alpha_s}{2 \pi} \sum_j \left( q_j \otimes P_{g \leftarrow q} +  g \otimes P_{g \leftarrow g} \right),$$ where $P_{\bar q \leftarrow q} = 0$ at leading order and $P_{\bar q \leftarrow g} = P_{q \leftarrow g}$ with $$(X_i \otimes P_{i \leftarrow j})(x) \equiv \int_{x}^1 \frac{d y}{y} X_i(y) P_{i \leftarrow j} \left(\frac{x}{y} \right).$$ Their addition gives $$\sum_i \mu^2 \frac{d}{ d \mu^2}  \left( q_i(x) + \bar{q}_i(x) + g(x) \right) = \frac{\alpha_s}{2 \pi} \left( \sum_i q_i \otimes (P_{q \leftarrow q} + P_{g \leftarrow q}) + g \otimes ( P_{g \leftarrow g} + 2 n_f P_{q \leftarrow g} ) \right)$$ and so $$ \frac{\alpha_s}{2 \pi} \int_0^1 \text{d} x \, x \, \left( \sum_i q_i \otimes (P_{q \leftarrow q} + P_{g \leftarrow q}) + g \otimes ( P_{g \leftarrow g} + 2 n_f P_{q \leftarrow g}) \right) = 0.$$ Now, noting that $$ \int_0^1 \text{d}x \int_x^1 \text{d}y = \int_0^1 \text{d}y \int_0^y \text{d}x,$$ and making a change of variables $x = yz$ decouples the $x$ and $y$ integrations allowing for the relations  given above for all $q_i$ and $g$. These relations are valid only at leading order, at the next order additional flavour effects come into play which makes the DGLAP evolution matrix structurally more involved. 
Additional comments:
The factor of $z$ included means it is the total momentum fraction that is conserved, not the number density. In this sense, it is $zP_{ij}$ that are the probability densities over $z$. Note the sum of probabilities of all possible splittings plus probability for no splitting is what sums to unity. You can trade the non integrable singularity at $z=1$ in the splitting functions for something integrable plus explicit IR divergences centred at $z=1$, this is done practically via the plus prescription. Schematically, this gives the following
$$\int_0^1 \text{d}z\, z \,( P^0_{g \leftarrow g} + P_{q \leftarrow g} + A_0 \delta (1-z) + \delta(1-z) ) = 1 \Rightarrow \int_0^1 \text{d}z\, z \,(P_{g \leftarrow g} + P_{q \leftarrow g} ) = 0, $$ where $P_{g \leftarrow g}$ is the plus-prescribed (regularised) splitting function and $A_0$ is of course tuned to cancel the divergence at $z=1$.
