Automated two-loop fermion-bubble contributions at high multiplicity

We compute the gauge-invariant, helicity-summed two-loop $N_f$ contribution for $q\bar q\to nγ^{(*)}$, $n=2,\ldots,8$. Local ultraviolet renormalisation, infrared factorisation, analytic loop-energy integration and a new automated threshold subtraction yield finite integrals over real spatial momenta. Adaptive multi-channel Monte Carlo samples loop momenta, helicities and identical-particle permutations. We quantify generation and evaluation costs across multiplicities. The same representation is integrated directly with parton distributions through five photons, yielding inclusive and differential virtual corrections.