Derivative-squared Yang-Mills action ($\mathcal{L}=\kappa (D_{\mu }{F}_{\nu \rho },D^{\mu }{F}^{\nu \rho })$)

This gauge-invariant higher-derivative density has equationThe extra commutator comes from varying the connection inside the adjoint covariant derivative: $\delta (D_{\mu }{F}_{\nu \rho })=D_{\mu }\delta F_{\nu \rho }+[\delta A_{\mu },F_{\nu \rho }]$. Two uses of gauge-covariant integration by parts give the first term, and Killing-form invariance gives the second. In an abelian limit it reduces to $\square {\partial }_{\mu }{F}^{\mu \nu }=0$. Keeping the covariant derivative order is essential because their commutator is curvature.