Factorization of squared amplitudes and the QCD splitting functions

When defining splitting functions from the Feynman rules, one makes use of the factorization of amplitudes. For definiteness, consider the process

$$a\,X \rightarrow b\,v\, X \rightarrow b\,W,$$

where parton $$a$$ splits into partons $$b$$ and $$v$$, and the latter collides with the target $$X$$, producing some set of particles $$W$$. In the limit where the virtual particle goes on-shell, the factorization is

$$\mathcal{M}_{aX\rightarrow bW} = \sum_{\rm \lambda_{v}} \mathcal{M}_{a\rightarrow bv} (\lambda_{v}) \frac{\rm i}{P_{v}^{2} - m_{v}^{2}} \mathcal{M}_{vX\rightarrow W}(\lambda_{v}) \,,$$

where $$\lambda_{v}$$ stands for the polarization of the virtual particle.

Next, the cross-section for the whole process is written as

$$d\sigma (a\,X\rightarrow b\,W) = \frac{1}{(\beta_{X} + 1) \, 2E_{a} \,2E_{X}} \,\, \overline{|\mathcal{M}_{aX\rightarrow b W}|^{2}}\, \frac{d^{3}{\bf p}_{b}}{2E_{b}} \, d\Phi_{W} \, \,,\\ \,=\, dz \,\, \frac{d{\bf p}^{2}_{\rm T}}{{\bf p}^{2}_{\rm T}} \,\, \frac{\alpha}{2\pi} \, P_{b\leftarrow a}(z) \,\, d\sigma (v\,X\rightarrow W)\,,$$

with the (unregularized) splitting functions are defined as

$$P_{b\leftarrow a}^{}(z) = \frac{1}{2 (2\pi)^{2}} \frac{z(1-z)}{p_{\perp}^{2}} \overline{|\mathcal{M}_{a\rightarrow bv}|^{2}}\,,$$

which means that the interference terms are disregarded in computing probabilities. Why is this done?

The first explanation that comes to mind is that the polarization of the virtual particle is (in the limit where it goes on-shell), definite but unknown, and thus one must sum over the squares of the amplitudes in computing the cross-section. However, it seems like this result should appear as a limit of the general case (where the intermediary particle is not close to being on shell).

The second possible solution that occurred to me is to simply compute the interference terms, which are proportional to $$\sum_{\rm \lambda_{a}, \lambda_{b}} \sum_{\rm \lambda_{v}} \mathcal{M}_{a\rightarrow bv} (\lambda_{v}) \cdot \mathcal{M}_{a\rightarrow bv}^{*} (-\lambda_{v})\,.$$

Note that, by the factorization property, and due to the fact that conjugating an amplitude turns incoming particles into outgoing particles (and vice-versa), this is proportional to

$$\sum_{\rm \lambda_{a}} \sum_{\rm \lambda_{v}} (P_{b}^{2} - m_{b}^{2}) \widetilde{\mathcal{M}}\Big(a(\lambda_{a}) v(\lambda_{v})\rightarrow a(\lambda_{a}) v(-\lambda_{v}) \Big)\,,$$

with the caveat that one of the vertex functions must be conjugated (which would only affect axial-vector couplings, I believe). The advantage of this result is that the sum is performed over configurations such that only one of the helicities is flipped. This means that, by angular momentum conservation, the amplitude must vanish for $$\theta=0$$, and thus must be at least of order $$\mathcal{O}(p_{\rm T})$$. Finally, since $$(P_{b}^{2} - m_{b}^{2}) \sim \mathcal{O}(p_{\rm T}^{2})$$, this produces an order $$p_{\rm T}$$ contribution to the splitting function, spoiling the large logarithm approximation.

Is any part of my reasoning valid? If not, is there a proper explanation of why the interference terms are neglected?