What are the gauge invariant parts of the Yang-Mills field strength tensor? In electromagnetism, which is an abelian gauge theory, we have the nice fact that all of the components of
$$
    F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu
$$
are gauge invariant quantities. We can equivalently talk about the field in terms of $\mathbf{A}_{\perp}$ and $\mathbf{E}_{||}$, the solenoidal part of the vector potential and the irrotational part of the electric field (respectively), as being a compete set of gauge invariant quantities that fully determines the state of the field everywhere without degeneracy.  
When we move to non-abelian Yang-Mills theory the field strength components
\begin{align}
    F_{\mu\nu} & = \frac{i}{g}\left[D_\mu,\, D_\nu\right] \\
        & = \partial_\mu A_\nu - \partial_\nu A_\mu - ig \left[A_\mu,\, A_\nu\right] \\
       & = T^a \left(\partial_\mu A^a_\nu - \partial_\nu A^a_\mu +gf^{abc} A^b_\mu A^c_\nu\right)
\end{align}
become matrices in the group space ($T^a$ are generators of transformations of the group). Thus the elements of the field strength indexed by space-time, $F_{\mu\nu}$, are no longer gauge invariant but gauge covariant — transforming as a rank two tensor in whatever representation of the gauge group the the generators are in.
When I want to know the group invariant parts of a matrix I think "eigenvalues", at least when the group is defined by preserving an inner product on some vector space (as $\operatorname{SU}(N)$ is). In other words, it seems to me like the physical information will be contained in the gauge invariant parts of $F_{\mu\nu}$ and finding the eigenvalues sounds like the way to go about it.
Has this been done? It seems unlikely that a single gauge transformation could simultaneously diagonalize all of the $F_{\mu\nu}$ at any given point in space-time, what are the consequences of this (if it's true)? Finally, would the eigenvalues of $F_{\mu\nu}$ also be independent of which representation the $T^a$ are in (i.e. would finding the eigenvalues in the defining representation give the same result as any other representation), assuming the representations all have the same normalization condition (e.g. $\operatorname{Tr}(T^a T^b)= \delta_{ab} / 2$)?
 A: The posed problem is formidable. If we had closed analytical expressions of the gauge invariant components of the reduced space of a gauge system, it would be a step towards the solution of the Yang-Mills Millennium problem. (We would then need to quantize this space, which is a formidable problem by itself).
Huebschmann, Rudolph and Schmidt  have performed the above exercise on a (very) simplified case of a lattice gauge theory consisting of a single plaquette; and as you can see, even this case is extremely difficult. The reduced phase space in this case is
$$M_{red} = T^*T/W$$
Where $T$ is the maximal torus of the gauge group $G$, $ T^*T$ is its cotangent bundle and $W$ is the Weyl group.
Moreover, the reduced space is not a manifold. It has the structure of what is called a stratified symplectic space, which is a disjoint union of strata of the action of the gauge groups on the unreduced space. The approximate computation of the energy levels in this problem entails the computation of tunneling probabilities between the strata.
Returning to the original question in the continuum
Some of the invariants of Yang-Mills fields are can be given by characteristic classes, please see the review by Eguchi, Gilkey and Hansen (section 6). These characteristic classes are gauge differential forms consisting of traces of polynomials of the Field strengths, which in addition, when integrated over appropriate cycles of the space-time manifold, give topological invariants. For example the second Chern class 
$$c_2 = \frac{1}{8 \pi^2} \left (\mathrm{tr} F \wedge F - \mathrm{tr}F \wedge \mathrm{tr}F \right )$$
However, these invariants do not separate the reduced gauge invariant space (i.e., do not exhaust the gauge invariant coordinates of the reduced space).
Another type of invariants based on holonomies involves the Yang-Mills connection (vector potential) and not the field strengths. These invariants have the form Wilson loops.
$$\mathrm{tr}\mathrm{P}e^{\int_{\Gamma}  A}$$ 
where the integration is over a closed path $\Gamma$ and $\mathrm{P}$ denotes path ordering. 
The collection of the Wilson loops over all possible paths in space-time does separate the reduced phase space. Actually, this representation is the basis of the loop formulation of the Yang-Mills theory, please see the following review by Loll. However, this approach poses great difficulties in quantization.
