Fermion field structure in non-abelian gauge theories I am trying to understand the structure of the fermions in non-abelian gauge theories. Disclaimer: my question might be very trivial (I suspect the answer could simply be "a change of basis"), but I would be grateful is someone could shed some light if there's something deeper lurking around the corner.
Let's consider the Yang-Mills Lagrangian 
$$ \mathcal{L} = -\frac{1}{4} (F_{\mu \nu}^a)^2 + \bar{\psi}(i \gamma^\mu D_\mu -m)\psi $$
where $ D_\mu = \mathbf{1} \partial_\mu - i g A_\mu^a t^a$, then $\psi$ needs to have both Dirac ($\mu$) and colour (a) degrees of freedom. I get confused when I change the regular approach to the problem and I am not sure if it is a real issue or if I'm just overcomplicating/overthinking this.
For instance, let's consider $SU(N)$ YM where the generators $t^a$ are $N \times N$ matrices, and there are $N^2-1$ of them. Therefore, the covariant derivative $D_\mu$ is an $N \times N$ matrix, and the index $a$ above runs over $N^2-1$ value. 
When we start by contracting $D_\mu$ with the $\gamma$-matrices, 
$$ (\gamma^\mu D_\mu)_{ij} = \delta_{ij} \gamma^\mu \partial_\mu - i g (\gamma^\mu A_\mu^a) (t^a)_{ij}  $$
one gets an $N \times N$ matrix of $4 \times 4$ matrices. The corresponding $4N$-component object this matrix acts upon is the $N$ Dirac spinors arranged in a column. 
However, notice that if we start in the following way instead:
$$ (D_\mu \psi)_i = \partial_\mu  \psi_i - i g A_\mu^a (t^a)_{ij} \psi_j $$
we get that the covariant derivative acts on $\Psi \equiv (\psi_1, \cdots, \psi_N)$. For all we know, $\Psi$ has no spinor structure since we haven't contracted with the gamma matrices yet. 
Contracting with the $\gamma^\mu$, one gets a $4 \times 4$ matrix of $N \times N$ matrices that encodes the same information as before. This time, it would seem that we only have one Dirac spinor, where each component is a $N$-valued singlet. 
However, it seems that the matrix $\gamma^\mu D_\mu$ and the $4N$-component spinor look different, even though all we did was to change the order in which we constructed things. 
In the first case, we get $N$ spinors corresponding to the $N$ colors of the adjoint representation. In $SU(3)$, this would be like saying we effectively have 3 spinors which correspond to the red/blue/green colors as $(\psi_R, \psi_B, \psi_G)$. In the second case, this identification fails since we only have one big complicated object. 
What went wrong? Is this difference simply a change of basis for $\psi$? Is there something relevant we can learn from looking at the YM Lagrangian in these two different ways? 
Also, in the first case, when we get $N$ spinors, I am confused about their signification. I always assumed that in QCD, $\psi$ would correspond to a quark, which is a fermion by itself. Does it mean that quarks are fermions that can be described by fermionic fields/degrees of freedom that we call color? 
 A: In a sense, you are right. It's the same information, put in a different order. Try to do that with a simpler case, for instance two two-dimensional vector space, one with index $\alpha=1,2$, the other with index $a=1,2$. Then any vector of this now four-dimensional vector space can be written either $\psi_{\alpha,a}$ or $\psi'_{a,\alpha}$ or even $v_i$ with $i=1,2,3,4$. Then the matrix $\sigma$ and $S$ acting respectively on the space $\alpha$ and $a$ can associated to create
$\gamma=\sigma \otimes S$ or $\gamma'=S\otimes \sigma$  acting respectively  on $\psi_{\alpha,a}$ and $\psi'_{a,\alpha}$, or any other complicated 4 by 4 matrix acting $v$ depending on how we associate a given $i=1,2,3,4$ to a set $(\alpha,a)=(1,1), (1,2),(2,1),(2,2)$. 
The choice to associate $i$ to $(\alpha,a)$ is free, but some choices are more natural than others. In you case, we are used to this of the system as $N$ fermions, described by spinors, which corresponds to $1\otimes \gamma_\mu$, because we think of the well defined spinors in interaction through the gauge field. But this is just a convenient representation.
A: You may write : 
$(\gamma^\mu D_\mu) (\psi_i)_k = [\delta_i^j  \partial_\mu - i g  A_\mu^a (t_a)_{i}^j] ~~ (\gamma^\mu)^k_l ~~(\psi_j)_l \tag{1}$
Here $i,j$ are in $1..N$, and $k,l$ are in $1..4$, $(\psi_i)_k$ is the k-th ($1 \leq k \leq 4$) component of the i-th ($1 \leq i \leq N$) spinor.
We could use the notation $\psi_{j~l} = (\psi_j)_l$, now we see that $[\delta_i^j  \partial_\mu - i g  A_\mu^a (t_a)_{i}^j]$ is acting on the first indice of $\psi_{j~l}$, while $(\gamma^\mu)^k_l$ is acting on the second indice of $\psi_{j~l}$, so we could use a tensorial compact notation, with $\psi$ representing the $\psi_{{j~l}}$ : 
$(\gamma^\mu D_\mu) ~\psi = ([\mathbb{Id}~ \partial_\mu- i g  A_\mu^a t_a] \otimes \tag{2}\gamma^\mu) ~\psi $
Of course, it is not very much useful for practical calculus, but it is, at least, a view on the whole structure.
