It seems to me the 'naive' approach to proving the existence of Yang-Mills in a rigorous context (via Osterwalder-Schrader $\to$ Wightman axioms), would be:
Study gauge invariant lattice QCD correlation functions such as $$\langle \mathrm{Tr}(W(\gamma_1)) \cdots \mathrm{Tr}(W(\gamma_n)) \rangle_{a,V} := \int dU e^{-S_\mathrm{Wilson}[U]} (\mathrm{Tr}(W(\gamma_1)) \cdots \mathrm{Tr}(W(\gamma_n))$$ $\gamma_i$ are physical paths in 4d, and $W(\gamma_i)$ are the corresponding wilson loops. The average is done on a lattice with lattice spacing $a$, and physical volume $V$. This observable is fully formal, finite, and something that I can go simulate on my computer right now. We hope that the 'Schwinger Functions': $$\langle \mathrm{Tr}(W(\gamma_1)) \cdots \mathrm{Tr}(W(\gamma_n)) \rangle := \lim_{a \to 0, V \to \infty} \langle \mathrm{Tr}(W(\gamma_1)) \cdots \mathrm{Tr}(W(\gamma_n)) \rangle_{a,V}$$ also exist (note that $\gamma_i$ are physical-sized loops, in other words they grow to contain more lattice links in the limit $a \to 0$)
Prove that the Schwinger Functions satisfy the Osterwalder Schrader axioms (and for the cherry on top, if the connected correlation functions have exponential decay then you have a mass gap)
What is the difficulty in this approach? For example, is it known whether the Schwinger Functions (step 1) even exist? It's certainly clear from numerical simulations that they should. If they are known to exist, step 2 doesn't seem too bad, given results like this.
(On physics.stackexchange, I've mostly seen discussion about other approaches, via algebras of observables and things like that. Sorry if this is a duplicate)