Solution (source code)

= Solution

The <axion> variation gives $\partial\mathcal L/\partial a=-m_a^2a+g_aF_{\mu\nu}\widetilde F^{\mu\nu}/4$, so
$$
\boxed{(\Box+m_a^2)a=\frac{g_a}{4}F_{\mu\nu}\widetilde F^{\mu\nu}.}
$$
For the <photon> variation, the <dual electromagnetic field tensor> obeys $\delta(F_{\rho\sigma}\widetilde F^{\rho\sigma})=2\widetilde F^{\rho\sigma}\delta F_{\rho\sigma}$. Antisymmetry then gives
$$
\frac{\partial\mathcal L}{\partial(\partial_\mu A_\nu)}=-F^{\mu\nu}+g_a a\widetilde F^{\mu\nu}.
$$
There is no undifferentiated $A_\nu$, hence
$$
\boxed{\partial_\mu F^{\mu\nu}=g_a\partial_\mu(a\widetilde F^{\mu\nu})
=g_a(\partial_\mu a)\widetilde F^{\mu\nu}.}
$$
The last equality uses the electromagnetic <Bianchi identity> $\partial_\mu\widetilde F^{\mu\nu}=0$. These signs follow from the positive $aF\widetilde F$ interaction in the source; reversing that interaction reverses both source terms. The induced electromagnetic current is identically conserved: its divergence contracts a symmetric second derivative of $a$ with an antisymmetric tensor, while the remaining term vanishes by the same <Bianchi identity>.