On the same spirit of this unanswered question I am proposing this question which I have been trying for some time now. Here I'm working with dimension $n = 4$ (identifying $\mathbb H = \mathbb R^4$) considering the principal $SU(2)$-bundle with $\rho : SU(2) \to GL(2, \mathbb C)$ being the adjoint representation $ad_g (A) = g^{-1} A g$, for all $g \in SU(2)$ and $A \in \mathfrak {su}(2)$. Moreover, the gauge potential is written (in local coordinates) as $\mathcal A = \mathcal A_{\alpha} dx^\alpha$, $\alpha = 1,2,3,4$ and its gauge field strength (curvature) is given by $$\mathcal F = d\mathcal A + \frac{1}{2}[\mathcal A, \mathcal A] = \frac{1}{2} \mathcal F_{\alpha\beta} dx^\alpha \wedge dx^\beta$$ where (after some lengthy calculation) $$\mathcal F_{\alpha\beta} = \partial_\alpha \mathcal A_\beta - \partial_\beta\mathcal A_\alpha + [\mathcal A_\alpha, \mathcal A_\beta]\ \ , \ \ \beta= 1,2,3,4$$ where $\partial_\alpha = \frac{\partial}{\partial x^\alpha}$.
Question:
Derive the Euler-Lagrange equations of the Yang-Mills functional given by $$\mathcal {YM} (\mathcal A) = \frac{1}{4}\int_{\mathbb R^4} \|\mathcal F\|^2 d(\bf vol_{\mathbb R^4})$$ such equations are called in the physics literature Yang-Mills equations $$\ast d^{\mathcal A}(\ast \mathcal F) = \sum_{\alpha = 1}^4 (\partial_\alpha \mathcal F_{\alpha\beta} + [\mathcal A_\alpha, \mathcal F_{\alpha\beta}]) = 0$$ where $\ast$ is the Hodge dual operator, $d^\mathcal A$ is the covariant exterior derivative.
Attempt: Since the space of gauge potentials is an affine space we may consider the variations $\mathcal A + t\mathcal B$ (family of 1-parameter of gauge potentials) thus
$$\begin{aligned}\frac{d}{dt} \Big(\mathcal {YM}(\mathcal A + t \mathcal B)\Big)\Big|_{t= 0} &= \frac{1}{4}\int_{\mathbb R^4} \frac{d}{dt} \Big(\|\mathcal F\|^2\Big)\Big|_{t = 0} d (\bf vol_{\mathbb R^4})\\ &=\frac{1}{2} \int_{\mathbb R^4} \|\mathcal F\| \frac{\left\langle \mathcal F, \frac{d}{dt} (\mathcal F)\Big|_{t =0}\right\rangle}{\|\mathcal F\|} d(\bf vol_{\mathbb R^4})\\&=\frac{1}{2} \int_{\mathbb R^4} \left\langle \mathcal F, \frac{d}{dt} (\mathcal F)\Big|_{t=0}\right\rangle d(\bf vol_{\mathbb R^4})\end{aligned}$$ where
$$\begin{aligned} \mathcal F &= \frac{1}{2} \mathcal F_{\alpha \beta} dx^\alpha \wedge dx^\beta \\&= \frac{1}{2}\Bigg( \partial_\alpha (\mathcal A_\beta + t \mathcal B_\beta) - \partial_\beta (\mathcal A_\alpha + t \mathcal B_\alpha) + [\mathcal A_\alpha, \mathcal A_\beta] \\&+ t [\mathcal A_\alpha, \mathcal B_\beta] + t [\mathcal B_\alpha, \mathcal A_\beta] + t^2 [\mathcal B_\alpha, \mathcal B_\beta]\Bigg) dx^\alpha \wedge dx^\beta\end{aligned}$$ by linearity we then obtain
$$\frac{d}{dt} (\mathcal F)\Big|_{t = 0} = \Bigg(\partial_\alpha \mathcal B_\beta - \partial_\beta \mathcal B_\alpha + [\mathcal A _\alpha, \mathcal B_\beta] + [\mathcal B_\alpha, \mathcal A_\beta]\Bigg) dx^\alpha \wedge dx^\beta$$ Now
$$\begin{aligned}\frac{d}{dt}\Big(\mathcal {YM}(\mathcal A + t \mathcal B)\Big)\Big|_{t = 0} &=\frac{1}{2} \int_{\mathbb R^4} \left\langle \mathcal F, \frac{d}{dt} (\mathcal F)\Big|_{t=0}\right\rangle d(\bf vol_{\mathbb R^4}) \\&= \frac{1}{2}\int_{\mathbb R^4} \mathcal F_{\alpha\beta} (\partial_\alpha \mathcal B_\beta - \partial_\beta \mathcal B_\alpha + [\mathcal A _\alpha, \mathcal B_\beta] + [\mathcal B_\alpha, \mathcal A_\beta])d(\bf vol_{\mathbb R^4}) \end{aligned}$$ Next step should be integration by parts. Then using the fact (Divergence Theorem)
$$\int_{\mathbb R^4} \partial_i (f) g dV_g = - \int_{\mathbb R^4} f \partial_i (g) dV_g $$ and that $\mathcal F_{\beta\alpha} = - \mathcal F_{\alpha\beta}$ we have
$$\begin{aligned}\delta\mathcal {YM}(A)&=\frac{1}{2} \int_{\mathbb R^4} \left\langle \mathcal F, \frac{d}{dt} (\mathcal F)\Big|_{t=0}\right\rangle d(\bf {vol}_{\mathbb R^4}) \\&= \frac{1}{2}\int_{\mathbb R^4} \mathcal F_{\alpha\beta} (\partial_\alpha \mathcal B_\beta - \partial_\beta \mathcal B_\alpha + [\mathcal A _\alpha, \mathcal B_\beta] + [\mathcal B_\alpha, \mathcal A_\beta])d(\bf {vol}_{\mathbb R^4})\\&= \frac{1}{2}\left(\int_{\mathbb R^4} \mathcal F_{\alpha\beta} (\partial_\alpha \mathcal B_\beta) d(\bf {vol}_{\mathbb R^4}) - \int_{\mathbb R^4}\mathcal F_{\alpha\beta}(\partial_\beta \mathcal B_\alpha) d(\bf {vol}_{\mathbb R^4})\right) \\&+ \frac{1}{2}\left(\int_{\mathbb R^4} \mathcal F_{\alpha\beta}[\mathcal A_\alpha, \mathcal B_\beta] + \mathcal F_{\alpha\beta}[\mathcal B_\alpha, \mathcal A_\beta])d(\bf {vol}_{\mathbb R^4}) \right)\\&= -\frac{1}{2}\left(\int_{\mathbb R^4} 2\mathcal B_\beta (\partial_\alpha \mathcal F_{\alpha\beta} d(\bf {vol}_{\mathbb R^4})\right) + \frac{1}{2}\left(\int_{\mathbb R^4} \mathcal F_{\alpha\beta}[\mathcal A_\alpha, \mathcal B_\beta] + \mathcal F_{\alpha\beta}[\mathcal B_\alpha, \mathcal A_\beta]d(\bf {vol}_{\mathbb R^4}) \right)\\&=-\int_{\mathbb R^4} \mathcal B_\beta (\partial_\alpha \mathcal F_{\alpha\beta} d(\bf {vol}_{\mathbb R^4}) \\&+ \frac{1}{2}\left(\int_{\mathbb R^4} \mathcal F_{\alpha\beta}[\mathcal A_\alpha, \mathcal B_\beta] + \mathcal F_{\alpha\beta}[\mathcal B_\alpha, \mathcal A_\beta])d(\bf {vol}_{\mathbb R^4}) \right)\\&=-\int_{\mathbb R^4} \mathcal B_\beta (\partial_\alpha \mathcal F_{\alpha\beta} d(\bf {vol}_{\mathbb R^4}) + \int_{\mathbb R^4} \mathcal F_{\alpha\beta}[\mathcal A_\alpha, \mathcal B_\beta] d(\bf {vol}_{\mathbb R^4})\end{aligned}$$
I need only show that
$$\int_{\mathbb R^4} \mathcal F_{\alpha\beta}[\mathcal A_\alpha, \mathcal B_\beta] d(\bf {vol}_{\mathbb R^4}) = - \int_{\mathbb R^4} \mathcal B_\beta [\mathcal A_\alpha,\mathcal F_{\alpha\beta}]d (\bf {vol}_{\mathbb R^4}) $$
Where do I want to get? $$\begin{aligned}\delta (\mathcal {YM}) (\mathcal A) = - \int_{\mathbb R^4} \mathcal B_\alpha (\partial_\alpha \mathcal F_{\alpha\beta} + [\mathcal A_\alpha, \mathcal F_{\alpha \beta} ]) d(\bf vol_{\mathbb R^4})\end{aligned}$$ then the stationary points satisfy
$$\sum_{\alpha = 1}^4 (\partial_\alpha \mathcal F_{\alpha\beta} + [\mathcal A_\alpha, \mathcal F_{\alpha\beta}]) = 0$$