Leibniz rule holds for covariant derivatives, both in gauge theories and gravity. Mathematically, a derivation is one for which the Leibniz rule holds. How does it work for non-abelian covariant derivatives. I will give you an example.
Let $\Phi^\dagger \Phi$ be invariant under local non-abelian gauge transformations. Then
$$
\partial_\mu (\Phi^\dagger \Phi) = (D_\mu\Phi)^\dagger \Phi + \Phi^\dagger (D_\mu\Phi) = [D_\mu(\Phi^\dagger)]\ \Phi + \Phi^\dagger (D_\mu\Phi)
$$
The term on the left hand side is the ordinary derivative as the combination is gauge invariant. By construction, it is easy to see that each term on the right hand side is invariant under local gauge transformation. You need to further show that the terms linear in the gauge field cancel. That more or less follows from representation theory. It is also easy to work out the rules for integrating by parts using the above formula. For your Lagrangian, you would consider the ordinary derivative of the scalar,i.e,, $\partial_\mu (\Phi^\dagger D^\mu \Phi)$.
Remarks:
- The ordinary derivative satisfies an additional property i.e., $\partial_\mu\partial_\nu\ f(x) = \partial_\nu \partial_\mu f(x)$ for all smooth functions. This is not assumed for all derivations. In fact, the obstruction to commutativity of the derivatives defines the field strength.
$$
[D_\mu,D_\nu]\ \Phi \sim g\ (F_{\mu\nu}^a T_a)\ \Phi\ .
$$
- In supersymmetry, one uses a graded version of the Leibniz rule for (covariant) derivatives involving Grassmann coordinates.
Let me first of all give a more precise definition of the curvature and exterior covariant derivative. To start with, lets fix the following data:
- A smooth manifold $\mathcal{M}$
- A principal $G$-bundle $P$, where $G$ is a (finite-dimensional) Lie group with Lie algebra $\mathfrak{g}$.
Now, as you have said, a connection $1$-form is a $\mathfrak{g}$-valued $1$-form, i.e. $A\in\Omega^{1}(P,\mathfrak{g})$ satisfying some extra properties. In order to define the corresponding curvature $2$-form $F^{A}\in\Omega^{2}(P,\mathfrak{g})$ there are several (equivalent) ways. One is to define it directly via the structure equation, which is
$$F^{A}:=\mathrm{d}A+\frac{1}{2}[A\wedge A]$$
where $\mathrm{d}$ is just the usual exterior covariant derivative of Lie-algebra valued forms, i.e.
$$\mathrm{d}A=\mathrm{d}(A^{a}T_{a}):=(\mathrm{d}A^{a})T_{a}$$
where $A^{a}\in\Omega^{1}(P)$ and $\{T_{a}\}_{a}$ is a basis of $\mathfrak{g}$ and where $[\cdot\wedge\cdot]$ is just the wedge product defined via the commutator, i.e.
$$[A\wedge A]:=\sum_{a,b}(A^{a}\wedge A^{b})[T_{a},T_{b}],$$
where $A^{a}\wedge A^{b}$ is just the standard wedge-product of real-valued forms. Starting from this, it is straight-forward to get the coordinate expression: Take a local section ("local gauge") $s\in\Gamma(U,P)$, where $U\subset\mathcal{M}$ open. Then, we define
$$A_{s}:=s^{\ast}A\in\Omega^{1}(U,\mathfrak{g})$$
as well as
$$F^{A}_{s}:=s^{\ast}F^{A}\in\Omega^{2}(U,\mathfrak{g}),$$
which are now forms defined on $U\subset\mathcal{M}$.
A straight-forward computation then yields
$$F_{\mu\nu}^{a}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A_{\mu}^{a}+c_{bd}^{a}A_{\mu}^{b}A_{\nu}^{d}$$
where $F_{\mu\nu}^{a}$ are defined via $F_{s}^{A}=F^{a}_{\mu\nu}T_{a}\mathrm{d}x^{\mu}\wedge\mathrm{d}x^{\nu}$ and $A_{\mu}^{a}$ via $A_{s}=A_{\mu}^{a}T_{a}\mathrm{d}x_{\mu}$. The constants $c^{a}_{bc}$ are the structure constants of the Lie algebra, i.e. $[T_{a},T_{b}]=c_{ab}^{d}T_{d}$.
Now, secondly, you can define the curvature via the covariant derivative $D_{A}$. Let $(\rho,V)$ be some (finite-dimensional) representation of $G$. Then, $D_{A}$ is a (family of) map(s) $$D_{A}:\Omega^{k}(P,V)\to\Omega^{k+1}(P,V)$$ defined via
$$(D_{A}\omega)_{p}(v_{1},\dots,v_{k}):=(\mathrm{d}\omega)_{p}(\mathrm{pr}(v_{1}),\dots,\mathrm{pr}(v_{k}))$$
for all $p\in P$, $v_{1},\dots,v_{k}\in T_{p}P$ and $\omega\in\Omega^{k}(P,V)$, where $\mathrm{d}$ on the right-hand side is the standard exterior derivative of $V$-valued forms (as explained above) and where $\mathrm{pr}:TP\to H$ is the projection onto the horizontal tangent space $H_{p}:=\mathrm{ker}(A_{p})\subset T_{p}P$ (the "Ehresmann connection corresponding to $A$").
Now, here is the most important point: You stated that $D_{A}$ can be
calculated via the formula
$$D_{A}\cdot=\mathrm{d}\cdot+\rho(A)\wedge\cdot.$$ That is in general not
true! In fact, it is only true for forms living in the subset
$$\Omega^{k}_{\mathrm{hor}}(P,V)^{\rho}\subset\Omega^{k}(P,V)$$
This set consists of all forms $\omega\in\Omega^{k}(P,V)$ satisfying the following two extra properties:
- "$\omega$" is horizontal": $\omega_{p}(v_{1},\dots,v_{k})=0$ whenever at least one of the tangent vectors $v_{i}$ is vertical (i.e. not contained in $H_{p}$).
- "$\omega$ is of type $\rho$": $(r_{g}^{\ast}\omega)_{p}=\rho(g^{-1})\circ\omega_{p}$ for all $g\in G$, where $r_{g}:P\to G,p\mapsto p\cdot g$ (the action $\cdot:P\times G\to P$ denotes the right-group action contained in the definition of a principal bundle)
In other words, your "definition" of the covariant derivative is actually only valid for this subset. Now, a connection $1$-form is in general not an element of $\Omega^{1}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$! (*)
To sum up, you cannot use your formula for a connection $1$-form. However, you can actually easily show that
$$D_{A}A\stackrel{!}{=}F^{A}=\mathrm{d}A+\frac{1}{2}[A\wedge A],$$
which then also leads to the correct coordinate expression as explained above.
* However, the curvature is an element of $\Omega^{2}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$ and also the difference of two connection $1$-forms is an element of $\Omega^{1}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$. The latter statement implies that the set of connection $1$-forms $\mathcal{C}(P)$ of a principal bundle is an infinite-dimensional affine space with vector space $\Omega^{1}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$.
Best Answer
You first need to know under which representation your field transforms, the flavor indices play no role here. Assuming $\phi$ is a scalar that transforms in the adjoint representation of the color group, sort of implied by the expression the OP writes. First we need to understand how the gauge covariant derivative acts on the field. Recall $T^a$ are infinitesimal generators of the Lie algebra, which means the only product that formally exists between them is the Lie bracket $$[T^a, T^b] = \sum_c {\rm i} f^{abc}T^c$$ where $f^{abc}$ are the structure constants. Then we have for example for a given flavor $f$, (which does not need to be manipulated) \begin{align} D_\mu \phi_f(x) &= \sum_a D_\mu^a \phi_f^a(x) T^a \\ &= \sum_a \partial_\mu \phi_f^a(x) T^a + \sum_{a,b} {\rm i}\,g\, \phi^a_f(x) A_\mu^b(x) [T^b,T^a]\\ &= \sum_a \partial_\mu \phi_f^a(x) T^a - \sum_{a,b,c}\,g f^{bac}\phi^a_f(x) A_\mu^b(x) T^c \end{align}
With that you should be able to work out the rest. The issue is, you need to identify the representation of your field. If $\phi$ transformed in the fundamental rep. then it would only have an extra index but no factors of $T^a$.