There are many similarities between the gauging of a $U(1)$ symmetry to obtain the physics of electromagnetism and the gauging of diffeomorphisms to obtain the physics of general relativity. In particular, in both cases one can obtain a covariant derivative: for electromagnetism,
$$ \nabla_\mu = \partial_\mu + i q A_\mu $$
where $A_\mu$ is the gauge field and for general relativity
$$ D_\mu = \partial_\mu + \Gamma_\mu^{\;ij}f_{ji} $$ where $\Gamma_\mu^{\;ij}$ is the connection and $f_{ji}$ is the generator of rotations.
However, these covariant derivatives behave quite differently: whilst in general relativity we have a nice 'distributive' property,
\begin{align} D_\mu a_\alpha b_\beta &= \partial_\mu a_\alpha b_\beta + \Gamma_{\mu\alpha}^{\;\;\gamma}a_\gamma b_\beta + \Gamma_{\mu\beta}^{\;\;\gamma}a_\alpha b_\gamma \\ &= (D_\mu a_\alpha) b_\beta + a_\alpha (D_\mu b_\beta) \end{align} which echoes the behaviour of normal derivatives, we don't have a similar relationship in electromagnetism:
\begin{align} \nabla_\mu a_\alpha b_\beta &= \partial_\mu a_\alpha b_\beta + i q A_\mu a_\alpha b_\beta \\ &= (\nabla_\mu a_\alpha) b_\beta + a_\alpha \partial_\mu b_\beta \\ &= a_\alpha (\nabla_\mu b_\beta ) + (\partial_\mu a_\alpha) b_\beta \\ &\neq (\nabla_\mu a_\alpha) b_\beta + a_\alpha (\nabla_\mu b_\beta) \end{align}
This makes me curious about the definition of the electromagnetism covariant derivative as a derivative at all. Does anyone have any insight they could share into the differences between these covariant derivatives, and in particular this lack of distributive property in electromagnetism?
I could imagine defining an operator $\mathcal{I}$ that has the same distributive property as $\partial_\mu$ and $f_{ji}$, and so $\nabla_\mu = \partial_\mu + i q A_\mu \mathcal{I}$ is distributive - is there a reason the covariant derivative isn't defined that way (or is it implicitly defined like that, and no one bothers to write it out fully)?