A gauge field transforms in the adjoint of the gauge group, but not in the adjoint (or any other) representation of the group of gauge transformations.
In detail:
Let $G$ be the gauge group, and $\mathcal{G} = \{g : \mathcal{M} \to G \vert g \text{ smooth}\}$ the group of all gauge transformations.
A gauge field $A$ is a connection form on a $G$-principal bundle over the spacetime $\mathcal{M}$, which transforms as
$$ A \mapsto g^{-1}Ag + g^{-1}\mathrm{d}g$$
for any smooth $g : \mathcal{M} \to G$. If $g$ is constant, i.e. not only an element of $\mathcal{G}$, but of $G$ itself, this obviously reduces to the adjoint action, so $A$ does transform in the adjoint of $G$, but not in the adjoint of $\mathcal{G}$. With respect to $\mathcal{G}$, it does not transform in any proper linear (or projective) representation in the usual sense, but like an element of a Jet bundle.
Let me discuss the mathematical precise derivation of the Yang-Mills equation in full generality. For this, let us first of all fix some notation:
- Let $P$ be a principal bundle over some (compact, oriented) pseudo-Riemannian manifold $(\mathcal{M},g)$ (spacetime) with structure group given by a compact (and finite-dimensional) Lie group $G$. The corresponding Lie algebra is denoted by $\mathfrak{g}$.
- Furthermore, let $\langle\cdot,\cdot\rangle_{\mathfrak{g}}$ be an $\mathrm{Ad}$-invariant inner product on $\mathfrak{g}$, or more generally, a non-degenerate $\mathrm{Ad}$-invariant and symmetric bilinaer form. For example, if $G$ is simple, then this is usually nothing else then a (negative multiple) of the Killing form.
Now, in order to formulate the action of Yang-Mills theory, we first of all take a connection $1$-form $A\in\Omega^{1}(P,\mathfrak{g})$, which corresponds to a gauge field in physics terminology. The corresponding curvature $F^{A}\in\Omega^{2}(P,\mathfrak{g})$ is defined by
$$F^{A}=\mathrm{d}A+\frac{1}{2}[A\wedge A].$$
Since we want to define an integral over the spacetime $\mathcal{M}$ and not $P$, we need to translate the curvature to a field defined on $\mathcal{M}$. This can be done in the following way: It is a general mathematical fact that there is the following isomorphism:
$$\Omega^{k}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}\cong\Omega^{k}(\mathcal{M},\mathrm{Ad}(P)),$$
where $\Omega^{k}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$ is some subset of $\Omega^{k}(P,\mathfrak{g})$ of forms satisfying extra properties (which are fulfilled by $F^{A}$) and where $\mathrm{Ad}(P)=P\times_{\mathrm{Ad}}\mathfrak{g}$ denotes the "adjoint bundle", which is a vector bundle defined on some certain quotient of $P\times\mathfrak{g}$. (In mathematics, this is a particular case of so-called "associated vector bundles", which can be defined for every principal bundle and every representation $(\rho,V)$ on $G$). As a second ingredient, we have to define an inner product
$$\langle\cdot,\cdot\rangle_{\mathrm{Ad}(P)}:\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))\times\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))\to C^{\infty}(\mathcal{M}).$$
This can be done in the obvious way: Take $\alpha,\beta\in\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$. If we take some local frame $\{e_{a}\}_{a=1}^{\mathrm{dim}(G)}\subset\Gamma(U,\mathrm{Ad}(P))$, which is a family of local sections of $\mathrm{Ad}(P)$, such that they form at each point a basis, then we can write
$$\alpha\vert_{U}=\sum_{a=1}^{\mathrm{dim}(G)}\alpha^{a}\otimes e_{a}\hspace{1cm}\text{and}\hspace{1cm}\beta\vert_{U}=\sum_{a=1}^{\mathrm{dim}(G)}\beta^{a}\otimes e_{a},$$
where $\alpha^{a},\beta^{a}\in\Omega^{k}(U)$ are real-valued forms. Using this, we define the iner product by
$$\langle\alpha,\beta\rangle_{\mathrm{Ad}(P)}\vert_{U}:=\sum_{a,b=1}^{\mathrm{dim}(G)}\langle\alpha^{a},\beta^{b}\rangle\langle e_{a},e_{b}\rangle_{\mathfrak{g}},$$
where $\langle\alpha^{a},\beta^{b}\rangle$ denotes just the usual inner product of real-valued forms and where $\langle e_{a},e_{b}\rangle_{\mathfrak{g}}$ has to be understood point-wise. After all this preliminaries, we define the Yang-Mills action via
$$\mathcal{S}_{\mathrm{YM}}[A]:=\int_{\mathcal{M}}\Vert F^{A}_{\mathcal{M}}\Vert^{2}_{\mathrm{Ad}(P)}\,\mathrm{d}\mathrm{vol}_{g}$$
where $F_{\mathcal{M}}^{A}$ is the curvature viewed as an element of $\Omega^{2}(\mathcal{M},\mathrm{Ad}(P))$ via the above isomorphism and where $\mathrm{d}\mathrm{vol}_{g}$ denotes the usual measure on pseudo-Riemannian manifolds. In order to derive the equations of motion, we first of all observe that
$$F^{A+t\alpha}=F^{A}+t(\mathrm{d}\alpha+[A\wedge\alpha])+\frac{t^{2}}{2}[\alpha\wedge\alpha]$$
for all $\alpha\in\Omega^{1}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}(P)}$. As a consequence, we have that
$$F^{A+t\alpha}_{\mathcal{M}}=F_{\mathcal{M}}^{A}+t(\mathrm{d}_{A}\alpha_{\mathcal{M}})+\mathcal{O}(t^{2}),$$ where $\alpha_{\mathcal{M}}\in\Omega^{1}(\mathcal{M},\mathrm{Ad}(P))$ corresponds to the form $\alpha$ via the isomorphism explained above. As a last ingredient, we need the following fact:
$$\int_{\mathcal{M}}\langle\mathrm{d}_{A}\alpha,\beta\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}=\int_{\mathcal{M}}\langle\alpha,\mathrm{d}_{A}^{\ast}\beta\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}$$
for all $\alpha\in\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$ and for all $\beta\in\Omega^{k+1}(\mathcal{M},\mathrm{Ad}(P))$, where $$\mathrm{d}^{\ast}_{A}:\Omega^{k+1}(\mathcal{M},\mathrm{Ad}(P))\to\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$$ denotes the codifferential. This formula is basically an extension of the well-knwon theorem for real-valued forms, that states that the exterior derivative is formally self-adjoint to the codifferential with respect to a suitable $L^{2}$-inner product, to bundle-valued differential forms. Using this, we get that
$$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{S}_{\mathrm{YM}}[A+t\alpha]=2\int_{\mathcal{M}}\,\langle\mathrm{d}_{A}\alpha_{M},F_{\mathcal{M}}^{A}\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}=2\int_{\mathcal{M}}\,\langle\alpha_{M},\mathrm{d}_{A}^{\ast}F_{\mathcal{M}}^{A}\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}\stackrel{!}{=}0.$$
Hence, by non-degeneracy of the inner product, which extends to the inner product $\langle\cdot,\cdot\rangle_{\mathrm{Ad}(P)}$, we get that
$$\mathrm{d}^{\ast}_{A}F_{\mathcal{M}}^{A}=0.$$
Using the definition of the codifferential, this is equivalent to say that
$$\mathrm{d}_{A}\ast F_{\mathcal{M}}^{A}=0,$$
which are the well-known Yang-Mills equation.
As a short comment: It is a general fact that for every connection $1$-form $A$ it holds that
$$\mathrm{d}_{A}F_{\mathcal{M}}^{A}=0.$$
This is a particular form of the "Bianchi identity". This equation together with the Yang-Mills equation can be viewed as a generalization of the well-known Maxwell's equations for electrodynamics.
Best Answer
Yes, classical plane wave solutions do exist for non-Abelian gauge theory. They were studied by Sidney Coleman; my example is equivalent to one of his under a gauge transformation. They are sometimes used as a "classical background" in research on the strong-field regime of QCD (but I don't know how the issue of quark pair production plays out). Coleman discusses plane waves, i.e., those with a single direction of propagation, whose field strength is constant in the orthogonal directions. But I see no problem with building a finite wavepacket, as long as we keep everything parallel to a single adjoint-representation vector $b_\alpha$.
As for whether gauge bosons are "charged", it's true that in these plane-wave solutions the currents, as classical observables, are zero. Their quantum expectation values also vanish (at least for one-particle states). But the current operators do have nonzero matrix elements between different one-particle states of this form. [EDIT: Of course, all talk of matrix elements for non-gauge-invariant operators such as the currents requires some gauge-fixing procedure. I think that a physical gauge choice, such as axial gauge $A^3_\alpha=0$, is simplest for my purpose here.]
So the most we can say, in the quantum theory, is that we have found a basis of one-particle states in which the currents have zero expectation. The existence of such a basis should not be surprising; we could always create a basis with this property, if desired, by taking superpositions of states with opposite charges. [EDIT: Actually, the fact that the expectation values for all of the $N^2-1$ charges vanish in the same basis does seem fairly surprising. My point was just that using this to call the bosons "uncharged" would be highly unjustified.]
The statement that gauge bosons are charged is fully justified by looking at each charge operator in its diagonal basis (that is, diagonal when restricted to the one-particle subspace). The matrices $A^\mu_\alpha$ in those basis states will be of rank two or more, meaning they will have different (probably orthogonal) adjoint-representation vectors for the different spatial directions, giving a nonzero current. In particular, the W bosons, which directly flip the weak isospin of fermions, do obviously have a nonzero value for one of the $SU(2)$ charges.