  We give a complete, self-contained, and mathematically rigorous proof that
Euclidean Yang-Mills theories are perturbatively renormalisable, in the sense
that all correlation functions of arbitrary composite local operators fulfil
suitable Ward identities. Our proof treats rigorously both all ultraviolet and
infrared problems of the theory and provides, in the end, detailed analytical
bounds on the correlation functions of an arbitrary number of composite local
operators. These bounds are formulated in terms of certain weighted spanning
trees extending between the insertion points of these operators. Our proofs are
obtained within the framework of the Wilson-Wegner-Polchinski-Wetterich
renormalisation group flow equations, combined with estimation techniques based
on tree structures. Compared with previous mathematical treatments of massless
theories without local gauge invariance [R. Guida and Ch. Kopper,
arXiv:1103.5692; J. Holland, S. Hollands, and Ch. Kopper, Commun. Math. Phys.
342 (2016) 385] our constructions require several technical advances; in
particular, we need to fully control the BRST invariance of our correlation
functions.
