# Non-Abelian gauge theory action with external currents: gauge invariance?

I am currently trying to work out the gauge invariance of the Yang-Mills action coupled to an external source and I noticed something weird.

The Yang-Mills action in presence of an external source $J^a_\mu$ reads

$S = \intop_x - \frac{1}{4} F^{a\mu\nu}F^a_{\mu\nu} + J^{a\mu} A^a_\mu = S_{YM} + S_{J}$.

The pure Yang-Mills part $S_{YM}$ is invariant under gauge transformations $V(x) \in \text{SU(N)}$ via

$\mathcal{A}_\mu = A^a_\mu t^a \rightarrow V \mathcal{A}_\mu V^\dagger - i V \partial_\mu V^\dagger$.

However, $S_J$, the part where we coupling to the external source happens, is not gauge-invariant. I know that the source should transform gauge covariantly via

$\mathcal{J_\mu} = J^a_\mu t^a \rightarrow V \mathcal{J}_\mu V^\dagger$,

but this does not render $S_J$ gauge-invariant. In fact I obtain

$S_J \rightarrow S_J - 2i\intop_x tr(V \mathcal{J_\mu} \partial^\mu V^\dagger)$,

which does not vanish. What's exactly happening here?

Note: if we look at the Abelian case and repeat the same steps we could actually make the action invariant by requiring charge conservation $\partial_\mu J^\mu=0$.

I know that the QCD action is gauge-invariant and in that case I obtain an extra term from the matter fields which cancels against terms from $S_J$. Does this mean that Yang-Mills with external sources is meaningless? Does Yang-Mills only make sense if the current $J^a_\mu$ comes from "dynamical" fields like quarks?

It is not possible to couple a conventional $c$-number source to a quantized non-abelian gauge field and maintain gauge invariance. For a current $J^\mu=\lambda^a J^\mu_a$, gauge invariance requires covariant current conservation $$0= \nabla_\mu J^\mu\equiv \partial_\mu J^\mu + [A_\mu,J^\mu]=0$$ and this requires a $c$-number(the first term) to equal a term containing an the operator $A_\mu$. It is because of this problem that external sources are introduced into a gauge theory as Wilson Loops $$W= P\exp\{ i\int A_\mu dx^\mu\}$$ where $P$ denoters a path-ordered integral of the matrix valued ($A_\mu = A_\mu^a \lambda^a$) gauge field.